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Abstract 

This  dissertation  documents  three  new  contributions  to  simulation  and  modeling 
of  optical  turbulence.  The  first  contribution  is  the  formalization,  optimization,  and 
validation  of  a  modeling  technique  called  successively  conditioned  rendering  (SCR). 
The  SCR  technique  is  empirically  validated  by  comparing  the  statistical  error  of 
random  phase  screens  generated  with  the  technique.  The  second  contribution  is  the 
derivation  of  the  covariance  delineation  theorem,  which  provides  theoretical  bounds 
on  the  error  associated  with  SCR..  It  is  shown  empirically  that  the  theoretical  bound 
may  be  used  to  predict  relative  algorithm  performance.  Therefore,  the  covariance 
delineation  theorem  is  a  powerful  tool  for  optimizing  SCR.  algorithms.  For  the  third 
contribution,  we  introduce  a  new  method  for  passively  estimating  optical  turbulence 
parameters,  and  demonstrate  the  method  using  experimental  data.  The  technique 
was  demonstrated  experimentally,  using  a  100  m  horizontal  path  at  1.25  m  above 
sun-heated  tarmac  on  a  clear  afternoon.  For  this  experiment,  we  estimated  « 
6.01  •  10-i)  m"i,  l0  17.9  mm,  and  L0  rs  15.5  in. 
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Modeling,  Simulation,  and  Estimation  of  Optical 

Turbulence 


1.  Introduction 

Applications  dealing  with  optical  turbulence  are  numerous  in  industry  and  the 
Department  of  Defense.  This  research  introduces  three  new  contributions  to 
help  advance  the  state  of  the  art.  The  following  sections  provide  a  roadmap  for  this 
document,  as  well  as  justification  for  research  it  represents  in  context  of  Air  Force 
needs. 

1 . 1  Organization 

The  following  is  a  list  of  chapter  descriptions: 

•  Chapter  1:  Serves  as  a  roadmap  for  the  document,  and  provides  justification 
for  this  research. 

•  Chapter  2:  Presents  a  modeling  technique  called  SCR  and  derives  a  theorem  to 
support  the  technique.  Empirical  data  are  used  to  validate  the  findings.  The 
contents  of  this  chapter  have  been  published  in  a  peer-reviewed  journal.  [9] 

•  Chapter  3:  Presents  a  method  for  estimating  optical  turbulence  parameters 
from  imaging  data.  The  method  is  demonstrated  using  physical  data  from  field 
experiments.  The  contents  of  this  chapter  have  been  submitted  to  the  journal, 
Applied  Optics,  and  are  currently  undergoing  peer  review. 

•  Chapter  4:  Provides  a  final  overview  of  this  research  by  reviewing  results  and 
conclusions  from  Chapters  2  and  3.  Outlines  the  direction  of  future  work  by 
this  author. 

•  Appendix  A:  Justifies  use  of  reetlinear  phase  integration  through  turbulent 
media  for  geometries  where  scintillation  (receiver  irradiance  fluctuation)  effects 
are  negligible. 
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1.2 


Contributions 


There  are  three  contributions  presented  in  Chapters  2  and  3: 

•  Contribution  1:  Formalization,  optimization,  and  validation  of  SCR  as  a  tech¬ 
nique  for  modeling  spatially  homogenous  and  isotropic  random  fields  in  n-dimensions. 
SCR  is  shown  to  be  a  statistical  filter  process  with  O(n)  computational  com¬ 
plexity,  where  n  represents  the  number  of  spatial  points  at  which  the  field  is  to 

be  simulated.  (Found  in  Chapter  2). 

•  Contribution  2:  Proof  of  the  covariance  delineation  theorem,  which  makes  pro¬ 
visions  for  computing  a  least  upper  bound  on  the  error  incurred  from  the  first 
level  of  SCR.  Its  use  precludes  any  need  for  costly  empirical  evaluation  while 
optimizing  a  particular  SCR  algorithm.  (Found  in  Chapter  2). 

•  Contribution  3:  Derivation  and  implementation  of  an  estimation  technique  for 
obtaining  optical  turbulence  model  parameters.  Finite  outer  scale  and  non-zero 
inner  scale  are  estimated  simultaneous  to  the  refractive  index  structure  constant. 

The  estimation  technique  is  demonstrated  using  physical  data  collected  from  a 
field  experiment.  (Found  in  Chapter  3). 

1.3  Value  to  the  Air  Force 

Contributions  1  and  2  are  good  for  rapidly  simulating  extensible  phase  screens 
to  an  arbitrary  level  of  detail.  Immediate  Air  Force  applications  include  characteriz¬ 
ing  performance  of  optical  sensors  and  algorithms.  Contribution  3  is  an  economical 
solution  for  characterizing  optical  turbulence  in  support  of  field  tests  of  Air  Force 
optical  systems.  Key  advantages  are  support  of  a  wider  range  of  propagation  path 
lengths  at  a  fraction  of  the  cost  for  scintillometer-based  alternatives.  Because  the 
sensing  technique  is  passive,  it  is  inherently  eye-safe. 
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1-4  Summary  of  Prior  Contributions  and  Theory 

Before  delving  into  the  details  for  the  new  contributions  presented  in  this  work, 
the  reader  may  wish  to  skim  this  section  as  a  primer  on  theory  related  to  previous 
works  published.  Each  subsection  discusses  material  as  it  relating  to  a  single  con¬ 
tribution.  Use  of  SI  units  is  implied  everywhere  in  this  document,  unless  specified 
otherwise. 

U,.l  Contribution  1.  This  contribution  consists  of  a  technique  for  simulat¬ 
ing  phase  screens  by  successively  conditioned  rendering.  Some  early  works  published 
on  phase  screen  simulation  relied  upon  Cholesky  decomposition.  The  recipe  for  this 
technique,  based  on  is  available  in  References  [24]  and  [25].  While  the  results  are  statis¬ 
tically  error  free,  the  Cholesky-based  approach  becomes  computationally  intractable 
for  applications  requiring  phase  screens  with  a  large  number  of  samples. 

In  1992,  Lane,  et  al,  used  a  midpoint  displacement  algorithm  from  the  computer 
graphics  community  and  applied  it  to  Kolmogorov  phase  screen  generation  [16] .  Lane 
showed  that  the  algorithm  was  significantly  faster  than  Fourier-based  techniques,  but 
less  accurate.  The  midpoint  displacement  algorithm  is  a  form  of  statistical  interpola¬ 
tion,  whereby  each  new  interpolate  is  obtained  by  adding  a  zero-mean  random  variable 
to  the  corresponding  interpolate  obtained  using  linear  interpolation.  In  this  paper, 
Lane  also  introduced  a  more  efficient  Fourier  technique  that  uses  a  non-uniform  sam¬ 
pling  grid  to  incorporate  subharmonic  frequencies.  Fewer  samples  are  used  overall, 
resulting  in  a  faster  simulation  than  the  one  presented  in  [25],  [24],  yet  similar  accuracy 
is  obtained  by  adding  more  samples  at  lower  spatial  frequencies.  The  downside  is  that 
the  fast  Fourier  transform  (FFT)  algorithm  cannot  be  used  with  nonuniform  sampling 
grids-a  fact  which  cuts  down  on  the  overall  computational  efficiency  of  the  method. 
As  a  followup  study,  in  1999,  Harding,  et  al  [13]  provided  a  more  in-depth  theoretical 
explanation  of  the  midpoint  displacement  algorithm  as  it  was  applied  to  the  problem 
of  phase  screen  generation.  Harding  also  examined  slight  variations  of  the  original 
algorithm,  and  characterized  the  computational  efficiency  of  the  algorithm.  In  2003, 
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Peterson  and  Mozurkewich  [23]  introduced  an  algorithm  that  applied  Lane’s  midpoint 
displacement  technique  to  selectively  interpolate  a  large-scale  phase  screen  for  mul¬ 
tiple  regions  of  interest.  They  also  added  the  concept  of  extending  the  phase  screen 
arbitrarily  in  any  direction.  By  contrast  to  prior  contributions,  the  work  presented 
in  this  document  offers  a  purely  covariance-based  Bayesian  probability  approach  to 
statistical  interpolation.  Furthermore,  a  variant,  of  Lane’s  algorithm  is  introduced, 
which  improves  both  the  efficiency  and  accuracy  for  this  class  of  algorithm. 

1J,.2  Contribution  2.  The  covariance  delineation  theorem  is  a  completely 
original  contribution  unrelated  to  prior  works.  Prior  to  introducing  the  theorem, 
we  provide  as  reference  a  few  definitions  specific  to  this  text.  These  definitions  are 
formally  introduced  in  Section  2.4: 

•  Correlated  Random  Set:  A  countable  set  of  correlated  random  variables. 

•  Conditional  Independence:  Statistical  independence  of  two  or  more  random  vari¬ 
ables  under  a  specific  set  of  conditions. 

•  Random  Estimate:  A  random  variable  whose  statistics  approximate  the  statistics 
of  another  random  variable. 

Readers  without  a  strong  background  in  linear  algebra  or  probability  and  statistics 
may  find  it  helpful  to  obtain  copies  of  References  [22]  and  [27]  to  facilitate  understand¬ 
ing  of  the  proof  for  this  theorem.  Both  texts  are  excellent  pedagogical  references  and 
use  similar  notation  as  this  document. 

1.4.3  Contribution  S.  This  contribution  is  a  new  approach  to  estimating 
optical  turbulence  parameters.  We  begin  by  defining  what  these  parameters  are  and 
discuss  their  historical  importance  in  optical  turbulence  modeling.  After  that,  we 
provide  a  synopsis  of  other  practical  methods  for  estimating  optical  turbulence  pa¬ 
rameters  found  in  the  literature. 
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1. 4-3.1  Optical  Turbule.nce  Parameters.  The  material  for  this  section 
comes  directly  from  Chapter  3  of  Laser  Beam  Propagation  Through  Random  Media 
by  Andrews  and  Phillips  [2].  Optical  turbulence  refers  to  spatial  fluctuations  in  the 
refractive  index  throughout  a  volume  of  atmosphere.  These  fluctuations  are  thought 
to  be  caused  primarily  by  temperature  fluctuations,  but  are  also  dependent  on  optical 
wavelength,  atmospheric  pressure,  and  relative  humidity  [4].  Such  fluctuations  occur 
at  random  and  with  some  amount  of  spatial  correlation.  We  refer  to  the  random 
spatial  variations  in  the  refractive  index  collectively  as  a  correlated  random  field.  A 
common  practice  in  modeling  optical  turbulence  is  to  assume  local  homogeneity  and 
isotropy  in  two  or  more  spatial  dimensions.  The  resultant  field  is  called  a  spatially  cor¬ 
related  homogenous  and  isotropic  random  field.  The  statistics  relating  points  within 
such  a  field  are  functions  only  of  the  distance  between  points,  and  implicitly  do  not 
depend  on  relative  orientation  or  absolute  position  of  points  within  the  field. 

In  1941,  Kolmogorov  developed  a  statistical  model  for  turbulence,  in  which 
he  related  the  statistical  structure  in  wind  velocity  measured  at  two  points  to  the 
distance  between  them.  The  turbulence  parameter,  L0,  was  used  to  denote  the  largest 
sized  eddy  with  homogenous  and  isotropic  statistical  structure.  Lq  is  known  as  the 
outer  scale  of  turbulence.  Another  parameter,  /,)•  was  defined  as  the  scale  for  which 
turbulent  eddies  lose  all  their  energy  to  viscous  forces,  lo  is  called  the  inner  scale  of 
turbulence.  Kolmogorov  defined  the  following  2/3  power  law  valid  for  distances  of 
separation  less  than  Lq  and  greater  than  Up 

Dv(R )  =E  [(V(P2)  -E[V(P2)})  -  (K(P.)  -E[V{P{) ]))2] 

-E[(V(P2)-1/(P,))2]  U-1) 

=c*rI, 


where 


R=\P-i-  Pi 


(1.2) 
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for  two  points,  P\  and  P2,  in  Euclidean  three-space.  Dy{R)  is  the  structure  function 
for  a  statistically  homogenous  and  isotropic  random  velocity  field,  E\-}  denotes  the 
ensemble  average  of  a  random  quantity,  V  represents  velocity,  and  Cy  is  the  velocity 
structure  constant. 

The  inertial  range  is  defined  as  a  subinterval  of  the  range  of  distances  between 
the  inner  and  outer  scales  of  turbulence: 

Zo  <  R  <  L0.  (1.3) 

Since  this  relationship  is  rather  ambiguous,  in  practice,  we  use  an  alternate  definition 
for  the  inertial  range  based  on  a  comparison  of  the  Kolmogorov  model  to  more  recent 
theory: 

r  ■  {ios>o  {^}  <  -1}  n  {iogi°  {^}  >  l5}  •  (L4) 

The  Kolmogorov  model  for  turbulence  is  only  considered  valid  for  the  inertial  range. 
In  1962,  Obukhov  applied  Kolmogorov’s  model  to  describe  the  statistical  structure 
of  temperature  fluctuations.  He  hypothesized  that  velocity  and  thermal  structure 
would  have  proportionate  statistical  structures,  because  the  energy  represented  by 
the  turbulent  flow  of  air  results  almost  entirely  from  fluctuations  in  temperature. 
External  fluctuations  in  pressure  can  also  induce  turbulent  air  flow.  For  purposes  of 
this  research,  wo  assume  an  environment  that  supplies  nearly  constant  atmospheric 
pressure  and  relative  humidity,  such  that  the  corresponding  effects  on  refractive  index 
fluctuations  are  negligible.  Obukhov  published  the  following  statistical  model  for 
temperature  fluctuations  on  the  inertial  subrange: 


Dt(R)  =  C (1.5) 

where  Dt(R)  is  the  structure  function  for  temperature  (T)  and  Cf-  is  the  thermal 
structure  constant.  Andrews  and  Phillips  use  the  following  simplified  model  to  relate 
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temperature  fluctuations  to  refractive  index  fluctuations: 


n(P)  =  1  +  79  xlO-6^,  (1.6) 

where  P(P)  and  T(P)  represent  pressure  and  temperature,  respectively,  at  some 
point,  P  in  Euclidean  tlnee-space.  A  wavelength  of  A  =  500  nm  is  assumed.  Relative 
humidity  is  not  specified  by  the  text  cited.  Units  of  pressure  for  this  expression  are 
millibars.  To  extend  Kolmogorov’s  2/3  power  law  to  refractive  index  fluctuations, 
one  has  to  make  additional  assumptions  and  approximations.  First,  we  express  the 
structure  function  for  the  refractive  index: 


Dn(R)  =£l[(n(P2)  —  n(P|))2] 

~c*m, 


(1.7) 


where  C'l  is  the  refractive  index  structure  constant.  Next,  we  apply  (1.7)  to  (1.6). 
Setting  K  =  79  x  10_6P,  we  have 


Dn(R)  =E 


=K2E 

=k2e 


p _ i_y 

\T(P2)  T(P)J 
:(t(pi}-t(p2))2-  ' 
(t(pot(p2))2  • 


(1.8) 


Define 


r(P)  =  r0  +  r1(P), 


(1.9) 


where  To  =  E[T\  is  the  mean  spatial  temperature  for  the  turbulent  volume  of  interest. 
Next,  we  assume  that  T  obeys  a  normal  distribution  with  E2[T ]  »  var[T],  and  express 
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approximate 


Dn(R)  «/v2 
=K2 
=K2 
—K2 


g[(r(P,)-T(P2)^ 

i?[(T(P,)T(P2))2] 

PAR) 

P[(T(P,)T(P2))2] 

_ PAR) _ 

P[((To  +  T1(P1))(To  +  T1(P2)))2] 

_ PAR) _ 

E  [: r( f  +  t?(Tx(Pi)  +  T,(P2))2  +  (r1(p,)T1(p2))2] 


(1.10) 


Since  the  thermal  fluctuations  are  expected  to  be  small  compared  to  the  mean  envi¬ 
ronmental  temperature  given  in  units  of  Kelvin,  we  can  write  |T0|  3>  |Tj(Pi)|,  and 
then  further  approximate 


PAR)  ~K 


=K 


■  PAR) 

p\n\ 
2  PAR) 

n 

=K2&1 

J0 


Therefore,  it  is  implied  that 


cl  =K^ 

in 


79  x  io-(,p 


Cl. 


(1.11) 


(1.12) 


The  covariance  spectrum,  <!>„(«),  for  the  refractive  index  in  a  statistically  ho¬ 
mogenous  and  isotropic  random  field  is  related  to  the  statistical  field  structure,  Dn(R), 

by 

oo 

Dn(R)  =  8tt  J  K2$„W  (l  -  "Ajt)  •  (L13) 

0 
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whore  k  denotes  scalar  spatial  frequency.  The  associated  spectrum  for  the  Kolmogorov 
model  within  the  inertial  range  is  given  by 

$„(«)=  0.033C£k"T.  (1.14) 

Several  years  later,  a  more  accurate  model  for  the  spectrum  of  refractive  index  fluc¬ 
tuations  emerged,  called  the  von  Karman  spectrum: 

*»(*)  =  0.033C*exp{-^|}  (k2  +  (2tt)2L0-2)-^  ,  (1.15) 

which  is  valid  for  0  <  R,  n  <  oo  and  incorporates  parameters  /0  and  L0  in  addition 
to  C%.  In  this  contribution,  we  use  derive  an  covariance-based  optical  technique  to 
estimate  all  three  parameters  of  this  model  simultaneously. 


1.4-3. 2  Alternative  Estimation  Methods.  The  majority  of  pre-existing 
turbulence  parameter  estimation  techniques  rely  on  measurements  of  scintillation. 
Scintillation  is  often  defined  as  the  fluctuation  in  received  irradiance  due  to  atmo¬ 
spheric  turbulence.  Scintillation  is  affected  by  all  three  parameters  (Iq,  Lq,  and 
C2)  [30].  The  scintillation  index  is  defined  as  [1] 


2  A  E[P] 
'  E2[I] 


-1, 


(1.16) 


where  /  represents  irradiance.  For  an  unbounded  plane  wave  propagating  through 
weak  turbulence  conditions  (a]  <  1),  we  have  [1] 


a)  .  «  1.23C^At:“, 

1  weak  n  7 


(1.17) 


which  is  also  known  as  the  Rytov  variance.  In  this  expression,  At:  is  the  propagation 
distance,  and  k  =  For  Rytov  variances  much  greater  than  unity,  (TjveHk  overesti¬ 
mates  af,  and  turbulence  fluctuations  are  said  to  be  saturated  [1].  Typical  methods 
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for  determining  Cft  and  Ln  use  this  <rjwcuk  to  describe  irradiance  fluctuations  received 
from  a  collimated  laser  propagating  through  weak  turbulence. 

As  turbulence  conditions  become  stronger  (o'juiealc  >  3),  the  inner  scale  1ms  an 
increasing  effect  on  the  scintillation  index.  The  most  recent  contribution  published 
uses  a  model  for  the  scintillation  index  that  is  valid  for  all  turbulence  conditions  for 
large-aperture  receivers.  The  model  takes  the  form  [30] 

L{hC'n-  D)  —  exp  |,7f„{/,ur9C}(^o,  Lu,  C~,  D)  T  aL{is,tuin}(^h  C*,  D) j  -  1,  (1.18) 

where  a\ fn{rlu,.yr}  is  the  variance  in  the  log  of  the  irradiance  due  to  large-scale  fluctu¬ 
ations,  and  ofn{/„u((}  describes  the  same  quantity,  except  for  small-scale  irradiance 
fluctuations.  D  is  the  diameter  of  the  aperture.  By  using  three  different  sized  aper¬ 
tures,  Vetelino,  et  al,  were  able  to  conduct  an  experiment,  in  which  they  obtained 
the  first  simultaneous  measurements  of  1$,  Lo,  and  Cft  [30].  In  order  to  obtain  enough 
optical  scintillation  to  measure  all  three  parameters,  a  1.5  km  path  was  used.  Three 
separate  receiving  apertures,  and  150  mw  of  1.55  /im  collimated  laser  radiation  were 
required  for  this  active  experiment.  Until  now,  there  has  been  no  other  attempt  to 
estimate  all  three  parameters  simultaneously. 

By  contrast,  this  document  introduces  a  method  of  simultaneously  estimating 
the  same  parameters  passively,  using  a  single  receiving  aperture,  requiring  an  optical 
path  length  of  only  100  m  through  turbulence.  Instead  of  using  scintillation  phe¬ 
nomenology,  as  does  the  majority  of  prior  literature,  this  technique  relies  on  theory 
derived  from  the  tilt  covariance  expressions  relating  a  one-dimensional  nonlinearly- 
distributed  array  of  point  sources. 


10 


2.  Optimized  Phase  Screen  Modeling  for  Optical  Turbulence 


Atmospheric  turbulence  is  a  fundamental  factor  contributing  to  the  imaging  qual¬ 
ity  obtainable  by  ground-based  astronomy.  Imaging  through  atmospheric  tur¬ 
bulence  is  a  process  commonly  modeled  by  inserting  a  random  phase  screen  directly 
over  the  aperture  of  the  imaging  system,  with  zero  or  more  additional  phase  screens 
distributed  along  the  optical  path.  Therefore,  a  key  component  to  modeling  astro¬ 
nomical  observations  is  the  means  to  quickly  and  accurately  simulate  random  phase 
screens.  In  their  recent  text,  Roggemann  and  Welsh  outline  a  method  for  simulating 
random  phase  screens  directly  by  using  a  phase  covariance  function  [24],  This  direct 
covariance  method  suffers  practical  memory  limitations  scaling  with  Q(n 4),  where 
n  refers  to  the  number  of  samples  per  axis.  In  1992,  Lane  et  al  introduced  a  fast 
and  practical  method  for  statistically  interpolating  existing  phase  screens  to  a  finer 
resolution  [16].  His  algorithm  is  based  on  the  midpoint  displacement  technique  which 
is  borrowed  from  the  computer  graphics  community.  Six  years  later,  in  collabora¬ 
tion  with  Lane  and  in  association  with  the  University  of  Canterbury  in  New  Zealand, 
Harding  et  al  optimized  the  algorithm  for  speed  and  accuracy  with  respect  to  the  size 
of  the  interpolator  [13].  Harding  et  al  also  presented  a  more  detailed  derivation  of  the 
algorithm.  Both  efforts  focused  on  applying  the  midpoint  displacement  algorithm  to 
Kolmogorov  phase  screens. 

2. 1  Roadmap 

This  chapter  builds  on  prior  works  by  presenting  an  alternative  statistical  in¬ 
terpolation  algorithm  which  is  useful  for  modeling  phase  with  a  finite  outer  scale. 
More  importantly,  this  chapter  introduces  a  theorem  which  gives  a  theoretical  basis 
for  the  statistical  interpolator,  and  simultaneously  generates  a  criterion  for  optimizing 
interpolators  by  lion-empirical  means.  While  this  chapter  applies  the  theorem  only  to 
random  phase  screens,  the  theorem  is  also  valid  for  any  other  simulation  that  imple¬ 
ments  SCR.  First,  we  provide  a  mathematical  introduction  to  SCR.  After  the  theorem 
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is  stated,  it  is  proven,  and  later  empirically  validated  by  comparing  two  alternative 
statistical  interpolators. 

2.2  Successively  Conditioned  Rendering 

Let  W  be  an  ordered  set  of  correlated  scalar  random  variables,  and  let  V  be  a 
corresponding  ordered  set  of  scalar-  realizations.  In  particular,  if  we  are  given  a  set 
of  realizations,  V  =  {vx,v2,  v2, . . . ,  %},  for  M  corresponding  correlated  random  vari¬ 
ables,  W  =  {wi,W2,w3, . . .  we  wish  to  simulate  one  realization  of  the  random 

variable,  w0  <£.  W,  by  rendering  a  corresponding  value  t>0  £  V.  After  v0  has  been 
assigned,  we  would  set  wM+i  <—  w0,  vM+1  <—  v0,  and  M  *-  M  - 1-1.  That  is  to  say, 
we  would  union  {wq}  with  W,  before  considering  a  new  random  variable  wq.  This 
process  is  referred  to  as  successively  conditioned  rendering  (SCR).  In  this  chapter  we 
are  concerned  with  random  phase  as  it  relates  to  imaging  applications.  Imaging  pro¬ 
cesses  are  inherently  insensitive  to  piston  and  residual  spatial  phase  is  predominately 
modeled  by  the  multivariate  normal  distribution  [24].  Therefore,  we  will  restrict  our 
SCR  analysis  to  the  zero-mean  multivariate  normal  distribution.  Let  v  be  a  column 
vector  containing  all  of  the  elements  of  the  set  V  and  let  w  be  a  corresponding  vector 
containing  all  the  elements  of  W.  The  joint  probability  density  for  the  set  of  random 
variables,  W,  defined  for  a  particular  set  of  realizations,  V,  is  given  by 

Pw(v)  =  [(27t)a/ |C|]_5  exp  {—^vTC~iv}  ,  (2.1) 

where  C  represents  the  covariance  matrix  describing  W,  and  is  defined  by  C  = 
E[(w  —  E[w})(u>  -  E[w})r\  =  E[wwT\.  We  define  w()  and  v0  to  be  the  column 
vectors  wu  =  [vc0, w'1  ] 7  and  v0  =  [u0,vt]T.  Similarly,  we  define  the  joint  pdf 
describing  «,'o  |J  W  as 

Pwu(v0)  4  [(27r)Af+1|Co|]-*  exp  ,  (2-2) 
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where  C0  =  E[w0w^].  Through  Bayes’  rule,  the  conditional  pdf  of  w0.  given  w  is 
defined  as 


(  \  _  Ptun(va) 

P™!"'(  °)  “  pjv) 

=  [(27t)m+1|Co|]~^  exp  {-^Cq  ^q} 

{(2n)M\C\]-hxp{-±yC-'v} 

=  Ki(v)  exp  {—\vq  Cq1v0 }  , 


where  we  have  conveniently  defined 


Ki{v)  ± 


\C\ 


l>|C0|J 


exp  {\vTC  ‘u} 


(2-4) 


to  encapsulate  all  the  coefficients  that  do  not  depend  on  wq .  Next,  we  define  So  — 
Cq1  —  {s,j},  where  i,j  €  {0, 1, 2, ,  M}.  Then,  the  quantity  VqCq1Vo,  taken  from 
the  argument  of  the  exponential  from  (2.3)  expands  to 


M  M 

vl'C^lv0  =  X  X  VjSij 

i=0  j= 0 
M  M 

=  EEwd  (2-5) 

«=o  o 

M  M 

=t>oSo,0  +  X  V'V0Si,0  +  X  V0VjS0j  +  K2(v), 

i= 1  3= 1 

where  we  have  isolated  the  terms  containing  v0  and  relegated  the  remaining  terms 
into  Ko(v).  By  definition,  C0  is  symmetric  [27].  Therefore,  its  inverse,  So,  is  also 
symmetric.  This  property  allows  us  to  combine  cross  terms,  resulting  in  the  expression 

M 

v'oCo'v  o  =  l’oS0,0  +  2  X  ViVoSifl  +  K2(v).  (2.6) 

i=l 
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Completing  the  square  would  require  a  conversion  of  the  form 


car  4-  bx  +  c  = 


i 

w-x 


2a  2 


+  c  — 

4a 


(2.7) 


substituting  x  <—  Vq.  Consolidating  the  remaining  terms  into  a  new  function,  K^(v), 
which  does  not  depend  on  vq,  enables  us  to  compactly  express  the  completed  square 
as 


/  M  x-1 

«o C~ 1  vn  =  I  sfov o  +  Sq|  Y  v’sc0  +  Ko(v) 


i=  1 
M 


=s|,o  Vo  +  so,0  Y  v,s,'°  +  Kz(v)- 


i=  1 


(2.8) 


The  resulting  form  for  the  conditional  pdf  from  (2.3)  is 


Pw„|i»(wo)  =  Ka(v)  exp 


M 


-550,0  Wo  +  «0,0  Y/  ViSi'° 


i= 1 


(2.9) 


which  we  immediately  recognize  as  a  univariate  normal  density  with  the  following 
parameters  for  the  mean  and  standard  deviation: 


M 

/to  =  —  s0,0  Y'j  ViSi,0 

»=i 


(2.10) 


From  this  form,  it  is  apparent  that  the  conditional  mean  is  simply  a  weighted  sum 
of  the  conditional  data,  and  that  the  conditional  standard  deviation  is  related  to  the 
first  diagonal  element  of  the  inverse  of  the  covariance  matrix.  For  applications  where 
C  is  reused  many  times,  it  will  suffice  to  pretabulate  the  first  row  of  C~' ,  allowing 
these  cached  values  to  be  reused  as  linear  filter  taps  for  obtaining  the  parameters 
defined  in  (2.10).  Under  such  conditions,  this  problem  reduces  to  complexity  O(M). 
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2.3  Statistical  Interpolation 

References  [13,16]  both  describe  the  resolution  enhancement  as  a  two-step  recur¬ 
sive  process,  where  the  first  step  is  linear  interpolation,  and  the  second  step  is  random 
midpoint  displacement.  This  chapter  introduces  a  new  perspective  on  achieving  the 
same  result.  The  two  steps  are  combined  into  a  single  step,  called  statistical  inter¬ 
polation,  and  is  based  on  SCR..  We  begin  with  a  set  of  M  phase  samples  uniformly 
distributed  over  a  coarse-resolution  2-D  rectangular  grid.  It  is  presumed  that  phase  is 
statistically  homogenous  and  isotropic  over  the  local  area  defined  by  the  phase  screen. 
This  assumption  allows  us  to  express  phase  covariance  as  a  function  of  separation 
between  two  points  in  the  aperture  without  regard  to  absolute  position  or  relative 
orientation.  We  define  B^(R±)  to  be  the  phase  covariance  between  two  points  in  the 
aperture  separated  by  the  distance,  /?x-  The  task  of  statistical  interpolation  is  to 
use  our  knowledge  of  existing  phase  samples  to  simulate  additional  phase  values  at 
locations  between  the  coarse-resolution  grid  points.  A  mathematically  pure  solution 
is  to  apply  SCR  to  render  each  additional  point  in  succession,  with  an  increasingly 
large  set  of  prior  data.  Unfortunately,  this  is  as  computationally  complex  as  applying 
the  direct  method  to  generate  a  fine-resolution  phase  screen. 

The  mathematically  tractable  approximation  originally  proposed  by  Lane  is  to 
limit  the  pool  of  prior  data  to  a  subset  of  V,  based  on  our  knowledge  of  the  covariance 
function.  The  phase  covariance  functions  for  all  published  turbulence  spectra  are 
positive  and  monotonically  decreasing  with  respect  to  R±.  This  property  implies 
that  the  nearest  neighbors  of  a  point  are  most  highly  correlated  in  phase.  This  is 
why  Lane’s  method  relies  on  the  subset  of  prior  data  represented  by  the  m  nearest 
neighbors  when  rendering  a  new  point  on  a  finer  resolution  grid.  Whereas  Lane 
and  Harding,  et  al,  pursued  this  idea  with  midpoint  displacement,  in  this  chapter, 
information  limited  SCR  is  applied  instead.  Otherwise,  the  algorithms  are  similar. 
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2-4  Covariance  Delineation  Theorem 

What  has  been  lacking  up  until  now  is  theoretical  proof  of  why  these  algorithms 
produce  such  accurate  phase  simulations  using  incomplete  information.  Furthermore, 
it  is  desirable  to  quantify  how  well  one  interpolation  scheme  is  expected  to  perform 
as  compared  to  another.  If  such  quantification  can  be  achieved  by  a  simple  formu¬ 
lation,  rather  than  by  mass  simulation,  one  would  be  able  to  quickly  optimize  such 
an  interpolator.  The  following  theorem  embodies  the  vital  link  between  theory  and 
practice. 

Definition  2.1  (Correlated  Random  Set).  A  countable  set  of  correlated  random 
variables. 

Definition  2.2  (Conditional  Independence).  Statistical  independence  of  two  or 
more  random  variables  under  a  specific  set  of  conditions. 

Definition  2.3  (Random  Estimate).  A  random  variable  whose  statistics  approxi¬ 
mate  the  statistics  of  another  random  variable. 

Theorem  2.1  (Covariance  Delineation  Theorem).  Consider  a  real-valued  cor¬ 
related  random  set,  W  =  {uq ,  W2,W:u  ■■.'%},  for  which  all  statistical  moments  are 
known.  Consider  also,  a  particular  corresponding  realization,  V  =  {<q ,  v2,  v:i,  ...vM}, 
of  the  correlated  random  set,  W.  Then,  consider  a  particular  correlated  random 
subset,  W'  =  {w\,W2,wl3, ...,  C  W,  with  corresponding  realization,  V7,  = 

{v\,V2,vtA, ...,  v‘m_\}  C  V,  which  is  collectively  referred  to  as  information  source  i. 
We  wish  to  examine  an  additional  real-valued  random  variable,  u>a  =  wo  £  W,  for 
which  all  statistical  moments  are  known,  but  for  which  no  known  realization  exists 
corresponding  to  the  members  of  Y.  The  statistically  correct  way  to  Simula te  u>n  is 
by  random  selection  according  to  the.  conditional  probability  density  function  (pdf). 
Pwt,\w(vh)\'io) ,  where  w  is  a  vector  containing  all  the  elements  of  W.  In  certain  ap¬ 
plications  where  computational  efficiency  is  at  a  premium,  it  may  be  more  desirable 
to  introduce  an  acceptable  level  of  error  by  restricting  the  amount  prior  data  used  to 
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condition  the  pdf  of  w0.  Consider  the  case  where  w0  is  approximately  simulated  ac¬ 
cording  to  the  conditional  pdf,  (wa\w'),  where  w1  is  a  vector  containing  all  the 

elements  of  W.  The  random  variable,  ib0,  described  by  this  pdf  is  a  random  estimate 
of  wo-  We  define  the  second-order  statistical  error  of  this  random  estimate  (SERE) 
as  the  normalized  root  mean  squared  error  between  the  second-order  statistics  of  Wq 
and  the  second- order  statistics  of  Wq  with  respect  to  W: 


£  = 


M 


X 


E[(w0  -  £[w0])2]  Y2  E^Wi  ~ 

i=0 
M 

Y2  I  (E[woWi]  -  E[w0]E[wj])  -  (E[w'0Wi\  -  E[wL0]E[wi\) 

i=0 

M 

E[(w0  -  E[w0])2}  Y2  E[(wl  -  £[wi]): 

«=o 
M 

Y2,  I  E[w0Wi\  +  E[wi\{E[w'u)  -  £[w0])  -  E[rb'0w, 


(2.11) 


L«=o 


Define  the  C  to  be  the  covariance  matrix,  C  =  {qj},  where  Cij  =  E[wLtvjj]  — 
E[wf\E[w'j\,  for  all  combinations  of  integers  {i,j}  €  [0,m  —  1]  x  [0,  m  —  1],  Then,  the 
least  upper  bound  on  SERE  that  is  computable  only  from  knowledge  of  C  is  given  by 


£<2[1-TJ*, 


(2.12) 


where  Tt  is  the  maximum  value  for  which 


Tt< 


Em— 1  r^m-1 

i=  1  2^>'=1  xiCijXj 


(2.13) 


is  true  given  any  set  of  real  values,  .  T,  is  called  the  statistical  degree  of  infor¬ 
mation  (SDI),  and  is  always  on  the  range,  [0,1],  with  Y,  —  1  corresponding  to  the 
case  of  perfect  prior  information. 
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2.5  Proof  of  the  Covariance  Delineation  Theorem 

We  shall  extend  C  to  include  an  mth  row  and  column  which  represent  the 
covariances  between  and  w0  U  W\  Let  the  extended  covariance  matrix  be  denoted 
by  C'.  Since  each  w,  except  w(>  was  jointly  simulated  according  to  pw(W),  and  w0  was 
approximately  simulated  according  top,i,0|W‘(u'o|W'),  the  only  undetermined  elements 
of  C"  are  dQm  and  dm0.  Since  C  is  a  covariance  matrix,  it  is  necessarily  symmetric. 
Therefore,  c^m  =  dm  0.  Thus,  we  really  only  have  one  undetermined  quantity  within 
C'. 


2.5.1  Positive  Semi- Definiteness  of  All  Covariance  Matrices.  Next,  we  will 
show  that  every  covariance  matrix  is  positive  semi-definite.  Tliis  proof  comes  from  a 
simple  extension  of  the  work  found  in  [22]  and  [27].  Consider  the  arbitrary  covariance 
matrix  C,  as  defined  above.  By  definition,  C  is  positive  semi-definite  if  and  only  if 


xTCx  >  0,Va:. 


(2.14) 


Expanding  this  product  into  a  sum  of  scalar  multiplications,  we  have  that 


m— 1  m-1 

x’  Cx  =  x’E\(wi  -  EKDK  -  E[wj]) ]xj 

i=0  j= 0 

‘m-1  m-1 

=  £  EE - EH])(vj - Eiwj\)xi 

.  ,=o  j= o 

‘m-1  m— 1 

=  e  J2  -  EKD  H  xiK-  - 

i=0  j= 0 


=  E 


>  0. 


( m-1 


Y.«-E  KD 


,  i=0 


(2.15) 


Therefore,  every  covariance  matrix  is  positive  semi-definite. 
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2.5.2  Continuity  of  Eigenvalues  in  Perturbations  of  a  Real  Symmetric.  Matrix. 

It  is  necessary  to  establish  the  property  of  continuity  of  the  eigenvalues  from  C'  with 

respect  to  symmetric  perturbations  in  elements  P0m  and  c'm  0.  Define  the  perturbed 
covariance  matrix.  Q  A  C1  +  A,  where  A0,m  =  A„li0  are  the  only  non-zero  elements  of 
A.  We  wish  to  show  that  the  eigenvalues  of  Q  are  continuous  in  A0,m.  Since  C'  and  A 
are  symmetric  and  real,  then  Q  is  also  symmetric  and  real.  Therefore,  the  eigenvalues 
of  Q  are  all  real.  The  eigenvalues  are  defined  as  the  roots  of  the  polynomial  |Q  —  A/|. 

This  polynomial  expands  to  the  form  | Q  -  A/|  =  a0  +  a]A  +  t^A2  H - 1-  am+iAm+1. 

Expanding  the  explicit  formula  for  the  determinant  found  in  [27],  will  reveal  that  each 
of  the  polynomial  coefficients  can  each  be  expressed  as  a  quadratic  polynomial  of  the 
variable  Ao,m.  Therefore,  the  coefficients  {a,}  are  continuous  in  A0>m.  The  roots  of 
a  polynomial  are  continuous  in  the  coefficients  of  that  polynomial  [5].  Therefore,  the 
eigenvalues,  {A,},  are  necessarily  continuous  in  the  coefficients,  {«,}.  Since  each  a, 
is  continuous  in  A<)_,„  and  each  A*  is  continuous  in  every  a*,  then  by  the  transitive 
property  of  continuity,  each  A,  is  also  continuous  in  Ao,m. 

2.5.3  Existence  of  an  Upper  Bound  on  the  Range  of  Positive  Semi- Definiteness. 

There  are  two  cases  to  consider  when  determining  the  allowable  range  on  A0,m  for 

which  Q  will  be  positive  semi-definite.  In  the  first  case,  we  restrict  C  to  be  positive 
definite.  For  this  case  we  wish  to  prove  that  there  exists  some  interval,  [A^rn,  A^„J, 
containing  zero  for  which  A0,m  results  in  Q  being  positive  semi-definite.  In  the  second 
case,  we  restrict  C  so  that  it  is  positive  semi-definite,  but  not  positive  definite.  For 
this  case  we  wish  to  prove  that  A0,m  =  0  is  a  necessary  condition  for  Q  to  be  positive 
semi-definite. 

First,  we  consider  the  case  where  C'  is  positive  definite.  We  define  the  com¬ 
ponents  of  some  non-zero  length  vector,  x',  according  to  x'  =  [xo,x\,X2,  ■  ■  ■  ,xmY  . 
We  also  define  x  A  [x0,xi,x2, . . .  ,zm_i]r  to  be  a  subvector  containing  the  first  m 
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components  of  x' .  Based  these  definitions,  we  can  write  the  following: 


x'tQx'  =  x/T(C'  +  A)*' 

=  x"  C'x'  +  x"  Ax' 

m  m 

=  x’c'x'  +  yy  yy  x,- a  ijxj 

i= 0  3=0 

=  x"  C'x'  +  2x0xmA0,m. 


(2.16) 


We  can  eliminate  the  trivial  case  where  xmx0  =  0,  since  it  corresponds  to  A(/m  —  Aq  m  - 
which  implies  that  the  error  is  unbounded  a  physical  impossibility.  Setting  x"  Qx'  = 
and  solving  (2.16)  for  A0,m  yields 


x^C'x' 


2 x0xm 

1 

m  m 

EE 

i= 0  j=0 

2x0Xm 

1 

m— 1  m— 1 

,  m— 1 

x  \  ^  ,  / 

r2  cJ 

E  E  XiCi’ixi 

i= 0  j= 0 

2  x0xm 

2^  X'Ci,mXm 

x°Xm  i=0 

2xQX7n 

x'Cx 

,  m— 1 

2X(jXm 

2x0 

A 

\C 

1— I  1 
1 

xm^Tn,m 

-  —  —  7  ZC'C  —  - 

2x0xm  x0  ^  4  2x0 


(2.17) 


where  we  have  implicitly  defined  A  =  x1  Cx  for  an  arbitrary  choice  of  x,  subject  to 
the  previously  established  constraint  that  ||x||  >  0.  Before  proceeding  any  further, 
we  will  divide  this  first  case  into  two  distinct  subcases. 


Consider  the  special  subcase  where  C'  is  positive  definite.  For  this  subcase,  we 


know  that  |Ao,m|  — 


:'rCV 


2xo*,. 


>  0,  which  precludes  the  range  of  Ao,m  from  containing 


values  within  some  finite  interval  about  A0,m  =  0.  Therefore,  a  finite  upper  bound 
exists  on  the  range  of  perturbation,  A0  tn,  for  which  Q  is  positive  semi-definite,  given 
that  C'  is  positive  definite. 


oo 
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For  the  other  subcase  C'  is  not  positive  definite,  but  it  is  still  positive  semi- 
definite.  For  this  subcase,  we  will  examine  the  function  A0im(a;m)  under  constant 
x.  This  function  is  continuous  everywhere  except  at  xm  =  0,  due  to  the  previously 
established  physical  requirement  that  xmx0  ^  0.  The  first  derivative  of  A0,m(arm)  is 
given  by 


dAo i7n  A  dmm 

dxm  2x0x'2m  2x0  ’ 

which  has  exactly  two  real  roots  defined  by 


(2.18) 


3’m  —  i 


A 


LLm,mJ 


(2.19) 


corresponding  to  the  critical  points  of  A0,m(:r,n).  Since  C  is  positive  definite,  we 
have  that  A  >  0.  dmm  is  also  positive,  since  it  represents  the  variance  of  a  random 
variable.  Consider  the  trivial  situation  where  dmm  =  0.  This  deterministic  constraint 
would  require  the  last  row  and  column  of  C'  to  be  all  zeros  to  maintain  positive  semi- 
definiteness.  Therefore  we  would  also  have  to  constrain  A0,m  =  0  in  (2.16)  for  Q  to 
be  positive  semi-definite.  For  all  remaining  situations,  we  must  have  dm  m  >  0,  since 
in  general,  dmm  represents  a  nonnegative  quantity.  Then,  for  A  >  0  and  drn  m  >  0, 


xh  = 


1  l 


<  0  and 


the  function  A0,m(:Em)  has  two  real-valued  critical  points,  .x“  =  -  |^A- 

A  1  >0,  on  opposite  sides  of  the  origin.  Next,  we  examine  the  second 


derivative  of  A0,m(xm), 


d2AQi)n(xm)  -A 


dxi. 


xox3m 


(2.20) 


which  tells  us  that  when  x0  >  0,  xfn  corresponds  to  a  local  minimum,  while  xbm 
corresponds  to  a  local  maximum.  The  reverse  is  true  for  x0  <  0.  This  information 
and  the  first  line  of  (2.17)  tell  us  that  for  xo  >  0,  we  have 


Ao, m(%m)  —  A|) ,m(xm)  >  0,Va:m  <  0 

Ao,m(^m)  <  A(l.m(x'|'n)  <  0,V:rm  >  0, 


(2.21) 
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while  for  Xq  <  0,  we  have 


Ao ,rn(xm)  <  A0,„, {x"n )  <  0,Vxm  <  0 

A o,m(zTO)  >  A0,m(^n)  >  0,Va;m  >  0. 


(2.22) 


Then,  regardless  of  the  sign  on  x0,  there  exists  some  interval,  (A£m,  A[/m),  for  which 
no  combination  of  x'  and  A0,m  can  be  found  to  satisfy  x"  Qx‘  =  0,  given  the  subcase 
of  positive  definite  C  and  positive  semi-definite  C‘ .  Therefore,  Q  is  positive  semi- 
definite  on  some  interval,  [A^m,  A^m],  where  A£m  <  A(,'m. 

For  the  second  case,  we  have  restricted  C  to  be  positive  semi-definite,  but 
not  positive  definite.  This  means  that  the  least  eigenvalue  of  C  is  zero  [27].  This 
corresponds  to  the  limiting  case  where  the  minimum  value  of  A  — >  0.  The  sign  of  A  is 
invariant  to  the  scale  of  x.  To  see  this,  we  substitute  x  =  ax  and  express  the  sign  of 
A  as 


A  x7  Cx 
]A|  =\xtCx\ 

( ax)7C(a.x ) 
|(ax)TC'(acc)| 
_  a2x‘ Cx 
|a2xTCi:| 
x1  Cx 
~\xTCx[ 


(2.23) 


which  is  clearly  unaffected  by  the  sign  or  magnitude  of  a  (a  >  0,  since  we  require 
||x||  >  0).  Therefore,  we  can  analyze  the  sign  of  A  by  restricting  x  to  the  set  of  unit 
vectors  without  loss  of  generality.  If  we  define  x  =  e,  where  e  =  [eo ,  e\,  dj, . . . ,  em-i]7 
is  the  unit  eigenvector  corresponding  to  Ao,  the  least  eigenvalue  of  C,  then  A  =  A0. 
For  Ao  >  0,  the  width  of  the  interval  of  A0,m  for  which  Q  will  be  positive  semi-definite 
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is  given  by 


A+  =|A0,m«)  -  A0,m« 


e0 


[*04,ml* 


Co 


H'H 


(2.24) 


2[Ao4 


N 


The  limit  on  this  expression  as  A0  — >  0  is  A+  =  0.  It  has  already  been  established 
that  this  interval  must  contain  zero.  Therefore,  as  A0  — >  0 ,  C  will  no  longer  be 
positive  definite  and  the  only  value  of  Ao,m  for  which  Q  will  be  positive  semi-definite 
is  A0,m  =  0,  which  proves  the  second  case. 


2.5.J,  Least,  Upper  Bound  on  Range  of  Positive  Semi- Definiteness.  In  Sub¬ 
section  2.5.3,  we  derived  an  expression  for  A+,  an  upper  bound  on  the  width  of  the 
interval  for  which  Ao,m  will  result  in  a  positive  semi-definite  Q.  In  this  subsection, 
we  will  minimize  A+  with  respect  to  x'  to  produce  the  least  upper  bound.  From  this 
point  forward,  we  will  only  consider  the  case  of  positive  definite  C.  since  we  have 
shown  that  there  is  zero  tolerance  for  error  when  C  is  not  positive  definite.  (2.19) 
provided  us  with  two  optimal  choices  for  xm  which  were  shown  to  produce  the  small¬ 
est  possible  A+  for  any  given  choice  of  x.  Then,  the  remaining  task  is  to  minimize 
A+  with  respect  to  x. 
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We  start  by  expanding  the  expression  for  A+  from  (2.24)  and  collecting  on  :r0 , 
writing 


A+  = 


2[xrCxdm>m\^ 


i= 0  j= 0 


2 

C'm,m  (  Y  Y  XiX&j  +  2X°  Y  Xi°0J  + 


i=l  j= 1 


(2.25) 


N  Cm’m  (  ^  °  ^  ^  :'?Wj  "f"  XqCo.0 


where  we  have  reverted  back  to  using  x  in  place  of  e,  and  we  have  redefined  A  = 
Yl"=\ 1  E7=i‘  xixjci,j-  Next,  we  examine  the  function, 


/(®  o)  = 


A  +  2;r0  Ejl/  XjCoj  +  3-qCo,o 


(2.26) 


which  is  defined  such  that 


A1”  —  2[cm,mJ{x  o))- 


(2.27) 


and  /(.To)  >  0  for  positive  semi-definite  C .  The  first  derivative  of  f(x0), 


has  only  one  root,  given  by 


df  2A  2  ^ 

o  o  /  ^j^O,j  ? 


■K0  v^m-1 


U  v^m-1 

Ej=l  XiC0,j 

The  second  derivative  of  /(:r0)  is  given  by 

d2/  6A 


(2.28) 


(2.29) 


(2.30) 
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When  evaluated  at  the  critical  point  given  by  (2.29),  this  expression  becomes 


d2f  2  1 

—  =-  ^xiCoJ  >0, 


0  lxn=sR 


(2.31) 


which  says  that  aig  corresponds  to  a  global  minimum  whenever  1  xjc oj  >  0. 
A  quick  analysis  of  (2.26)  tells  us  that  for  xjco,j  =  then  f(x o)  =  A  +  c0)o- 

Either  way,  .r0  =  ref,  minimizes  f(x o).  Performing  this  substitution  yields  the  following 
expressions: 


'  OJ  |x0=x!! 


xjcoj  , 

'  -A  ~+C».0 

XjCo.j 


—  CQ,0 


^  "m- 1 
-  Y^xjeo j 
Li=i 


(2.32) 


i  r-1 

=  2  cm  m(°o,o  1  :  ’y  '  XjCo,j 

ACQ.O  ~T 

L;=i  J 


(2.33) 


From  (2.33),  it  is  apparent  that  the  maximum  allowable  range  of  perturbation,  A+, 


will  be  minimized  by  maximizing  expression, 


^  ’m— 1  '  2 

T(*)  =  T—  . 

Ac°'°  YU 


(2.34) 


where  we  have  defined  the  a:  to  be  a  subvector  of  x,  such  that  x  =  [aq ,  X2,  £3, . . . ,  xm^i]T . 
Y,.  is  taken  to  mean  the  maximum  value  of  T(x),  and  is  called  the  statistical  degree 
of  information  (SDI).  Recall  that  A  =  Y^j=\  xixjCij,  so  that  T(x)  expands  to 


T(x)  = 


{K  ■:*,«*? 

Era— 1  \—^m—  1 

1=1  2^j=l  XiXjCi  j 


(2.35) 
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Then,  maximizing  the  T(i),  will  simultaneously  produce  the  least  upper  bound  on 
the  range  of  positive  semi-definiteness  of  Q 


2.5.5  Existence  of  the  Least  Upper  Bound.  Next,  we  will  show  that  a  single 
global  maximum  exists  for  T(x)  with  respect  to  each  xk  in  x.  Rewriting  T  as  a 
function  of  a  particular  xk ,  gives  us  the  form 


T(*fc)  =c0io 


m— 1  m—  1  m—1 

xW,k  +  2xk  £(1  -  Sjik)xjCjik  +  -  hk)il 

j=l  x=  1  j= 1 

m—1  ^ 

xkCo,k  +  53  (i  -  h*)xjc0,: 


=  (Ax\ 

\Co,oCk,kJ 


J=1 

gtixkY 


'h.^Xi'X'jCjj 


(2.36) 


where 


9i(*k)  =  \Xk  +  —  -  Sj^XjCoj  ) 

\  Co'fc  i= i  J 


m—1 


m—1  m—1 


<Mxk)  =x\  +  2-^  51(1  -  Sjik)XjCjik  +  —  Y]  V(1  -  5j,k)(  1  - 

Cfc-fc  i=i  Cfc'fc  «=i  j=i 


(2.37) 


For  brevity,  we  can  rewrite  <j\  and  go  as 


9i(xk)  =  (xk  -  hi)2 
92{xk)  =  {xk  -  b2)2  +  d2, 


(2.38) 
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where  bi,  bo,  and  d2  are  given  by 


b\  =  ~ 
h  =  ~ 


m—  1 


1  '  ‘ 

°'k  j=  1 

771 —  1 

—  hk)xici,k 


^  m— 1 m- 1 

d2  =  -b\  +  —  ^2  X^1  ~  5J.*)(1  “  h^XiXjCij. 
G*-*  i=l  J=1 


(2.39) 


From  these  forms,  we  recognize  cji  and  g2  to  be  concave  up  parabolas,  each  having 
a  single  positive  global  minimum.  We  wish  to  characterize  T (xk),  which  is  directly 
proportional  to  the  ratio  of  gt  and  go  with  a  positive  constant  of  proportionality.  We 
will  consider  two  cases  in  analyzing  this  T (x*). 

For  the  first  case,  we  consider  6,  =  b2.  According  to  (2.39),  the  condition 
b]  =  bo  would  require  the  trivial  solution,  x  =  0,  or  else  for  Hn  integers 

j  G  [l,m  —  1).  The  trivial  solution  is  excluded  as  it  presents  as  a  singularity  in  T. 
The  other  solution  can  be  interpreted  as  forcing  all  columns  of,  C,  the  lower-right 
principal  submatrix  of  C,  to  be  linearly  dependent.  The  implication  is  that  every 
element  in  column  k  of  C  would  need  to  be  defined  as  c,-,*  =  c*,*.  Furthermore, 
since  C  is  a  covariance  matrix,  it  must  be  symmetric,  forcing  Cjj  =  ck  k.  This  would 
make  every  element  of  C  identical.  Because  of  the  dependency  between  {cj0}  and 
the  columns  of  C,  every  element  of  C  would  also  be  identical.  Therefore  C  would 
be  positive  semi-definite,  but  not  positive  definite.  From  (2.35),  we  see  that  this  case 
corresponds  to  T  =  1,  or  perfect  information.  Therefore,  T  is  maximized  for  any 
choice  x  whenever  b\  =  b2. 

For  the  second  case,  we  consider  the  imperfect  information  case,  b\  d  &2-  Since 
A  >  0  and  c0(q  >  0,  we  must  also  have  that  (P  >  0,  Therefore,  T(x*)  is  strictly 
nonnegative  with  no  singularities  or  other  discontinuities.  Furthermore,  an  upper 

C2 

bound  on  the  range  of  T(xk)  is  and  from  (2.36)  and  (2.38),  we  see  that  as 
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xk  ->  ±00,  T(xfc) 


k  ■  The  first,  derivative  of  xk  is  given  by 

dT  f  co,k  \  <j2fl\  -  giffz 

dxk  \coficktkJ  gl 

. (  C<U  \  f  (('■'' k  ~  1'2)2  +  (f)(2) (xk  -bx)  (xk  -  bi)2{2)(xk  -  b2)\  .  , 

\co,oCk,kJ  \  ((xk-b2y+<py  b2y  +  (py>  )  [ '  1 

_  (  2c0,fc  \  / 4(^1  -  62)  +  xk(b \  -  b'i  +  (p)  ,  1)11)2  -  Inbl  -  6 1  r/2  \ 

Vco.oCfc,  J  v  ((**  -  62)2  +  d2)2  T  ((**  -  62)2  +  d2)2/  ’ 

which  has  the  two  real  roots, 


xk  =[2(£>i  -  62)]  1  [fc'f  -  62  -  d2  ±  [(62  -b\-  d2)2  -  4  (&[  -  62)(6?&2  -  bxb\  -  M2)]* 

=[2(6!  -  6a)]-1  [b‘i  -  b\  -  d 2  ±  [bj  +  b\  +d4  +  %\b2  -  4{b%  +  6, +  6,M2) 

+  2(6262  +  6?d2  +  62d2)]*' 

62-6  |-d2±[(61-b2)a  +  d2] 

2(61-62) 

(2.41) 

Because  T(x/..)  is  continuous  and  bounded  on  (—00, 00)  and  asymptotically  approaches 
unity  at  either  extreme,  one  of  these  roots  corresponds  to  a  global  maximum,  while  the 
other  corresponds  to  a  global  minimum.  To  determine  which  critical  point  corresponds 
to  the  maximum,  we  consider  two  subcases.  In  the  first  subcase,  we  let  bx  >  62.  Then, 
we  have  bx  >  b2  —  We  will  examine  the  sign  of  the  slope  for  regions  outside  of 
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the  interval,  jb2  —  -  b  1  ■  For  =  61  +  e,  for  any  positive  e,  we  have 

=  (~  ■--■■■)  [(61  -  62) 2  +  d2]-2  x  [(26,e  +  e2)(^  -  b2)  +  e{b]  -  b]  +  d2)] 
\C0fiCk,k/ 

=  f-—k-)  [(&i  ~  h)2  +  d2]  ~2  x  [(62  -  2bib2  +  b\  +  d2)e  +  (&,  -  b2)e2] 

=  f  ■  *"  y  ')  [(&i  -  b2)2  +  d2]  x  [((&!  —  b2)2  +  d2)e  +  (fri  -  62)e2] 

\  ca,oCk,k  / 


5T 


Xk=b\+e 


>0, 


(2.42) 


which  tells  us  that  xk  =  61  will  minimize  T (xk),  and  consequently,  that  xk  =  b->  - 
will  maximize  T  for  >  b>.  Conversely,  when  61  <  b2,  we  have  that  <  b2  —  , 

so  we  examine  for  a:*  <  61.  With  similarity  to  (2.43),  for  xk  =  b\  —  e,  given  any 
small  positive  e,  we  have 


dr 


dxk 


xk=bi-t 


+  d2)e  —  (b2  —  di)e2 


(2.43) 


which  tells  us  that  xk  =  bi  will  minimize  T  (a*) ,  and  consequently,  that  xk  =  b 2  — 
will  maximize  T  for  61  <  b2.  Therefore,  regardless  of  which  is  the  greater  of  61  and 
62,  xk  —  b2  —  will  maximize  T(xk).  Since  a  single  computable  maximum  value 
exists  for  r(xk),  given  any  set  of  elements  {a-,}  from  x,  subject  to  i  ±  k,  then,  a 
global  maximum  value  also  exists  for  T(x)  for  the  set  of  11011-zero  length  sub  vectors, 
x.  Substituting  this  global  maximum  for  Y(x)  into  (2.34)  and  then  into  (2.33)  is 
therefore  guaranteed  to  produce  the  least  upper  bound,  A+in  <  A+: 


^mtn  —  2[co,oCm,m (1  ^t)]  2  • 


(2.44) 
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Proof  of  convergence  for  an  iterative  algorithm  to  find  T  is  beyond  the  scope  of  this 
document. 

2.5.6  Upper  Bound  on  the  Second-Order  Statistical  Error  of  p,vu\w,.(wo\w'). 
Without  explicitly  using  our  knowledge  of  {Um0,  dmX,cJm2, . . . ,  it  is  impossi¬ 

ble  to  determine  the  location  of  the  center  for  the  allowable  range  of  Ao,m.  Therefore, 
we  must  assume  the  worst  case  scenario  where  C'  is  not  positive  definite,  and  conse¬ 
quently,  A0,m  must  be  either  strictly  nonnegative  or  strictly  nonpositive.  In  this  case, 
the  maximum  absolute  error  in  covariance  for  element  d0m  will  be  A+in.  Rewrit¬ 
ing  the  maximum  absolute  error  in  covariance  between  w0  and  a  particular  variable, 
Wi  tf.  W\  in  terms  of  expected  values  yields 


Krun  =  m(n>0  -  £M)2]£[(^  -  £H)2](1  -  TJ]i  (2.45) 


Using  this  error  bound,  the  normalized  mean  squared  statistical  error  of  the 
random  estimate  (SERE),  is  given  by 


IA+  |2  I  * 

^  Z— /i=o  I  t, mm  I 

-  [E[(w0  -  EH)2]  ZZo  E[(m  -  EM)2]. 

_  'm«>0  -  £K])2]  E"0g[K  -  E[Wi ])2](1  -  T,) 

Eliwo-EiwomZZoEliWt-ElwW] 

=2[1  —  TJ5, 


(2.46) 


which  is  sufficient  to  prove  the  covariance  delineation  theorem. 


2. 6  Experimental  Validation  of  the  Covariance  Delineation  Theorem 

The  covariance  delineation  theorem  provides  the  least  upper  bound  on  the  statis¬ 
tical  error  associated  with  a  particular  interpolator  at  any  level  of  recursion.  However, 
since  the  theorem  assumes  we  have  perfect  knowledge  of  the  covariance  statistics  for 
all  random  variables  under  consideration,  for  the  time  being,  our  analysis  must  be 
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restricted  to  the  first  level  of  recursion.  Successive  levels  would  rely  on  imperfect 
realizations  of  random  variables,  for  which  the  exact  covariances  are  unknown.  We 
shall  now  consider  two  possible  statistical  interpolators.  The  first  is  Lane's  4-nearest- 
neighbors  interpolator,  which  is  depicted  in  Figure  2.1.  The  second  interpolator  is 


Figure  2.1:  4-Nearest-Neighbors  Statistical  Interpolator.  This  figure  illustrates  the 

t.liree  generations  rendered  during  each  level  of  recursion  by  the  4-nearest  neighbors  recur¬ 
sive  interpolator.  Curved  arrows  indicate  the  cycle  of  generations,  while  the  black  arrows 
indicate  which  prior  data  (tails)  is  used  to  render  additional  data  (arrowheads).  Step  0 
represents  the  initial  coarse  grid  of  phase  samples. 


called  the  2-nearest-neighbors  interpolator,  and  is  depicted  in  Figure  2.2.  Both  inter¬ 
polators  render  the  first  two  generations  without  corrupting  the  input  data  that  gets 
used  for  each  new  sample  that  is  rendered.  Both  interpolators  rely  on  statistically 
imperfect  data  from  prior  generation(s)  to  render  the  third  generation.  The  key  differ¬ 
ence  between  these  two  interpolators  is  the  first  generation.  The  4- nearest- neighbors 
method  uses  four  neighbors  each  at  a  relative  distance  of  1.41  units  from  the  sample 
being  rendered;  by  contrast,  the  2-nearest-neighbors  method  uses  two  neighbors  each 
located  1  unit  away  from  the  sample  being  rendered.  In  the  second  generation,  both 
interpolators  use  two  neighbors  at  1  unit  separation,  while  for  the  third  generation, 
both  interpolators  use  four  neighbors  at  1  unit  separation.  The  obvious  question  con- 
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Figure  2.2:  2-Nearest-Neiglibors  Statistical  Interpolator.  This  figure  illustrates  the 

three  generations  rendered  during  each  level  of  recursion  by  the  2-nearest  neighbors  recur¬ 
sive  interpolator,  using  the  same  graphical  conventions  as  Figure  2.1 


cerning  the  first  generation  is  whether  it  is  better  to  use  two  neighbors  within  close 
proximity,  or  four  neighbors  that  are  40%  more  distant.  (We  are  also  given  that  both 
interpolators  in  this  scenario  use  circularly  symmetric  inputs.)  To  answer  this  ques¬ 
tion,  the  first  experiment  was  simulated  at  various  magnitudes  of  scale  with  respect 
to  Lq,  the  outer  scale  of  the  atmosphere. 

While  it  is  a  worthwhile  exercise  to  validate  the  theorem  with  the  first  exper¬ 
iment  using  only  first-level  first  generation  statistics,  it  is  of  more  practical  value  to 
consider  an  experiment  involving  phase  screens  constructed  from  multiple  levels  of 
interpolation.  While  the  stated  theorem  does  not  directly  apply  to  such  a  case,  one 
would  intuitively  expect  that  an  interpolator  with  greater  accuracy  in  the  first  gener¬ 
ation  of  the  first  level  of  recursion  implies  greater  statistical  accuracy  across  all  levels 
and  generations.  This  is  the  focus  of  the  second  experiment. 

Both  experiments  were  based  on  the  von  Karman  spectrum  of  optical  fluctu¬ 
ations.  The  inner  scale  was  set  to  zero,  and  the  Rytov  approximation  for  a  plane 
wave  in  turbulence  was  applied  to  obtain  the  following  phase  covariance  function  for 
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a  single  layer  of  turbulence,  published  in  [24] : 

OO 

B^R±)  =  0.033(2ir)2k2C2Az  f  J?MkR^  dK,  (2.47) 

l  (*  +  %)* 

where  k  is  the  wave  number,  C2  is  the  mean  refractive  index  structure  constant  for 
the  path  of  turbulence,  and  Az  is  the  total  path  length  within  turbulence.  It  is 
further  assumed  that  C2  varies  slowly  in  the  direction  of  the  optical  axis.  The  results 
published  in  this  chapter  have  been  normalized  with  respect  to  L0,  C2,  Az,  and 
wavelength,  A. 

2. 7  Results 

2.7.1  Experiment  1:  Single-Level,  Single- Generation  Interpolation.  For  the 
first  experiment,  7.5  million  2x2  phase  screens  were  simulated  and  statistically  in¬ 
terpolated  to  level  1,  generation  A  for  each  algorithm.  For  the  4-nearest-neighbors 
method,  this  means  using  four  points  to  simulate  a  fourth  point  at  some  distance, 
b,  from  each  of  the  priors,  as  depicted  in  Figure  2.1.  For  the  case  of  the  2-nearest 
neighbors  method,  it  means  using  the  leftmost  two  points  of  a  2x2  phase  screen  to 
simulate  a  fourth  point  at  a  distance  of  from  the  two  priors  used.  See  Figure 
2.2  for  the  depiction.  SCR  is  used  for  both  algorithms;  the  midpoint  displacement 
algorithm  is  not  used  in  this  work. 

A  maximum  likelihood  estimator  was  derived  for  estimating  the  covariances 
between  the  first  interpolated  point  and  the  four  original  points.  Figure  2.3  represents 
the  results  from  this  experiment.  The  figure  clearly  indicates  that  the  2-nearest- 
neighbors  interpolator  performs  considerably  better  than  the  4-point  interpolator  for 
the  first  level  and  generation.  This  was  predicted  by  the  ordering  of  the  least  upper 
bounds  for  the  statistical  error  of  each  interpolator.  The  2-point  interpolator  had  a 
slightly  lower  upperbound,  so  we  would  expect  it  to  perform  better  than  the  4-point 
interpolator. 
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Figure  2.3:  Statistical  Error  Estimates  for  Level  1,  Generation  A.  7.5  Million  simula¬ 
tions  were  performed  for  each  interpolator  to  obtain  maximum-likelihood  estimates.  Solid 
regions  indicate  99.9%  confidence  intervals. 


2.7.2  Experiment  2:  Multilevel  Interpolation.  The  focus  of  the  second  ex¬ 
periment  was  to  test  whether  the  theorem  has  predictive  value  in  determining  which 
interpolator  will  perform  better  across  multiple  levels  of  recursion,  in  the  presence  of 
accumulative  statistical  error.  For  this  experiment,  4.6  million  33x33  phase  screens 
were  interpolated  from  2x2  phase  screens.  Figures  2. 4-2. 6  indicate  the  same  trend 
that  was  predicted  by  the  covariance  delineation  theorem.  This  initial  experiment 
indicates  that  the  theorem  ma.y  be  extended  as  a  tool  for  optimizing  interpolators 
across  multiple  levels  of  recursion. 

2. 8  Conclusions 

In  this  chapter,  a  new  approach  to  statistical  interpolation  was  introduced,  based 
on  limited-information  successively  conditioned  rendering  (SCR.)  theory.  The  result¬ 
ing  algorithm  combines  the  previously  published  two-step  algorithm  into  a  single  step 
algorithm.  A  new  theorem,  called  the  covariance  delineation  theorem,  was  introduced 
as  a  means  for  optimizing  the  statistical  accuracy  of  SCR-based  applications,  such  as 
statistical  interpolation.  The  theorem  was  rigorously  proven  to  work  for  the  first- 
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»  in'*  Statistical  Error  for  Multilevel  Interpolation 


Figure  2.4:  Statistical  Error  Attributed  to  Each  Level.  This  plot  indicates  the  degree 
of  statistical  error  per  level  of  statistical  interpolation.  The  data  represent  4.6  million 
phase  screen  simulations  per  interpolator.  The  original  phase  screens  were  2x2,  resulting 
in  33x33  phase  screens  after  5  levels  of  recursion.  Dotted  lines  indicate  99.9%  confidence 
intervals. 
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Figure  2.5:  Histogram  of  Multilevel  Statistical  Error.  This  histogram  compares  the 

distribution  of  statistical  error  for  each  interpolator,  without  regard  to  any  particular  level 
or  generation.  The  99.9%  confidence  interval  has  a  mean  width  of  ±1.56%  of  the  estimated 
SERE  values  shown.  This  histogram  reflects  the  same  data  shown  in  Figure  2.4. 


level,  first  generation  stage  of  statistical  interpolation.  Two  different  interpolators 
were  tested  empirically  against  the  theorem,  and  the  results  were  consistent  with  the 
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Maximum-LIkefihood  Multilevel  Covariance  Estimates 


Figure  2.6:  Multilevel  Phase  Covariance  Estimates.  This  plot  compares  the  ideal  the¬ 
oretical  phase  covariance  function  to  the  actual  covariance  resulting  from  each  interpolator. 
This  plot  reflects  the  same  data  shown  in  Figures  2.4  and  2.5. 


theory.  A  second  experiment  was  performed  to  test  the  plausible  use  of  the  theorem 
for  optimizing  statistical  interpolation  across  multiple  levels  of  recursion.  The  results 
promising  enough  to  warrant  further  theoretical  study  of  this  problem. 
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3.  A  New  Approach  to  Estimating  Optical  Turbulence 

Statistics 


The  effects  of  optical  turbulence  have  been  studied  since  the  early  1960’s.  [1] 
Andrews  and  Phillips,  et  al  have  published  a  brief  treatise  [l]  on  the  matter 
as  it  relates  to  the  phenomenon  of  scintillation,  which  is  complemented  by  a  more 
complete  pedagogical  discussion  in  their  classic  text,  Laser  Beam  Propagation  Through 
Random  Media  [2],  Nearly  all  published  works  on  the  topic  of  characterizing  optical 
turbulence  rely  on  three  particular  underlying  assumptions  which  are  all  absent  from 
the  derivation  presented  in  this  chapter.  The  first  assumption  is  that  either  strong  or 
weak  turbulence  conditions  apply,  requiring  two  separate  analyses  of  the  problem  with 
regard  to  scattering  phenomenology,  and  hence,  scintillation.  The  second  assumption 
is  that  the  Rytov  approximation  will  lead  to  sufficiently  accurate  results.  Finally, 
and  perhaps  most  importantly,  the  published  solutions  for  Rytov-based  scintillation 
theory  all  employ  the  paraxial  approximation  to  optical  wave  propagation.  [2]  While 
there  are  many  other  assumptions  common  to  prior  works,  this  work  addresses  and 
circumvents  each  of  the  aforementioned. 

3. 1  Roadmap 

This  chapter  begins  with  a  review  of  prior  contributions  found  in  the  literature. 
Next,  we  provide  a  mathematical  derivation  which  begins  with  a  model  for  refrac¬ 
tive  index  fluctuations  and  ends  with  a  computationally  tractable  expression  for  tilt 
covariance.  Section  3.4  gives  a  description  of  how  the  experiment  was  conducted, 
what  observable  data  were  used  and  how  the  data  were  processed.  The  section  also 
provides  a  description  of  the  environmental  conditions  at  the  time  of  the  experiment. 
Results-orientcd  readers  may  wish  to  skip  to  Sections  3.5  and  3.6  which  present  and 
discuss  the  results  from  the  experiment. 
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3.2  Background, 

Several  prior  works  have  boasted  capability  to  measure  optical  turbulence  statis¬ 
tics,  but  all  are  limited  to  the  three  previously  stated  assumptions.  Three  statisti¬ 
cal  parameters  traditionally  used  to  characterize  turbulence  are  the  refractive  index 
structure  constant  (C£),  the  outer  scale  of  turbulence  (L0),  and  the  inner  scale  of 
turbulence  (lu).  In  1973,  Dunphy  and  Kerr  published  temporal  statistics  of  for 
a  low-altitude  6  km  path,  using  Kolmogorov-Rytov  scintillation  theory.  [6]  In  1979, 
Valley  working  in  conjunction  with  the  US  Air  Force  published  his  work  on  the  effects 
of  a  non-zero  Z0  and  finite  L(]  on  short-term  and  long-term  Strehl  ratios.  [29]  In  1985, 
Ochs  and  Hill  [21]  published  early  measurements  of  l0,  using  a  dual-scintillometer 
aparatus  originally  proposed  by  Livingston.  That  work  was  also  partially  sponsored 
by  the  US  Air  Force.  Two  years  later,  Frehlich  published  the  conclusion  that  a  single 
(paraxial)  phase  screen  does  not  adequately  model  inner  scale  effects  on  scintillation 
observations.  [10]  In  1988,  Coulman,  et  al  published  a  method  for  estimating  L0  using 
C'n  profiles  “obtained  from  spatio-angular  correlation  measurements  of  stellar  scintil¬ 
lation.”  [7]  The  aparatus  used  is  for  making  those  observations  is  called  scintillation 
detection  and  ranging  (SCIDAR).  Hill  published  another  work  in  1988  comparing 
two  different  scintillation  methods  for  estimating  l0,  where  the  dependence  of  C*  is 
cancelled  out  by  taking  the  ratio  of  two  quantities  that  are  proportionate  to  C;t  and 
dependent  on  ln.  [14]  Particular  advantages  were  identified  for  each  method.  In  1989, 
Martin  and  Flatte  used  a  multi-phase  screen  technique  to  simulate  the  scintillation  of 
a  spherical  wave  emanating  from  within  the  turbulence,  as  a  means  to  validate  exper¬ 
imentally  observed  phenomena.  [17]  The  paraxial  approximation  and  Rytov  method 
were  employed  between  each  phase  screen,  and  deterministic  spatial  filters  were  ap¬ 
plied  at  each  screen  to  allow  phase  screens  of  finite  extent.  In  1992,  Hill  and  Ochs 
published  a  work  in  which  micrometerological  measurements  of  Z0  were  in  agreement 
with  theory,  and  scintillation-based  estimates  of  l0  also  agreed  fairly  well.  This  work 
was  sponsored  in  part  by  the  US  Army.  Hill  and  Frehlich  later  published  on  related 
work  under  sponsorship  of  NASA,  NSF,  and  the  US  Army,  in  which  they  used  a 
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spherical-coordinate  scintillation  simulation  to  characterize  the  relationship  between 
l()  and  irradiance  variance  under  strong  turbulence  conditions.  [15]  In  the  following 
year,  Masciadri  and  Vernin  proposed  an  angle  of  arrival  technique  for  estimating 
l0  which  essentially  uses  tilt  variance  estimates  made  from  simultaneous  observa¬ 
tions  through  multiple  apertures  of  varying  diameters.  [18]  In  2000,  Flatte  and  Ger¬ 
ber  applied  multi-phase  screen  spherical-wave  and  plane-wave  simulations  to  the  Hill 
spectrum  of  turbulence  and  concluded  that  under  strong  turbulence  conditions  their 
results  did  not  match  first  order  analytical  predictions  of  scintillation.  That  work  was 
supported  in  part  by  the  US  Navy.  [8]  In  2002,  Whiteley  presented  an  analytical  tech¬ 
nique  for  using  a  multi-aperture  or  multi-source  system  to  estimate  a  Cl  profile.  [32] 
While  the  technique  is  not  scintillation-based,  the  derivation  still  applies  the  paraxial 
assumption  and  uses  the  Kolmogorov  spectrum.  That  work  was  sponsored  by  the  US 
Air  Force.  In  2004,  Ziad  and  Scliock,  et.  al ,  compared  three  astronomical  techniques 
for  measuring  L0,  including  two  methods  based  on  the  Shack-Hartmann  wavefront 
sensor,  and  a  third  interferometric  technique.  [34]  Also  in  2004,  Cain  published  a 
technique  under  Air  Force  sponsorship  for  C'l  profile  estimation  using  a  wavefront 
sensing  array  and  a  single  laser  guide  star.  Most  recently,  in  2006,  Vetelino  et  al 
published  results  from  a  2005  experiment  in  which  simultaneous  estimates  of  C'l,  /0, 
and  L0  were  computed  from  scintillation  observations  collected  using  three  separate 
receiving  apertures  of  different  sizes.  [30]  That  work  was  also  partly  supported  by  the 
US  Navy. 

In  this  work,  we  present  an  alternative  approach  to  the  problem  of  estimating 
statistical  turbulence  parameters.  Rytov  theory  will  not  be  applied,  and  neither  will 
the  paraxial  (a.k.a.  phase  screen)  approximation.  The  technique  being  explored  is 
unaffected  by  scintillation,  and  consequently  the  effects  of  atmospheric  scattering  can 
be  ignored.  The  technique  is  similar  to  the  angle  of  arrival  (AA)  methods  described 
fry  Masciadri  [18]  and  Ziad  [34],  and  to  the  differential  tilt-based  technique  proposed 
by  Whiteley  [32],  but  it  is  perhaps  most  closely  related  to  the  method  proposed  by 
Cain  [3],  and  is  in  actuality  an  extension  of  the  same  idea.  Instead  of  extracting 
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a  C'l  profile  from  tilt  variances  estimated  over  multiple  apertures,  we  are  using  a 
single  aperture  to  estimate  tilt  covariance  for  multiple  pairs  incoherent  point  sources. 
The  estimated  tilt  covariances  which  are  subsequently  applied  to  estimate  the  spatial 
covariance  of  refractive  index  fluctuations.  Since  the  resulting  system  of  equations 
is  underdetermined,  we  constrain  the  problem  by  applying  a  three-parameter  von 
Karman  model  with  finite  outer  scale  and  non-zero  inner  scale.  For  this  work,  we  chose 
a  one-dimensional  model  for  the  spectrum  of  refractive  index  fluctuations,  defined 
under  the  typical  assumptions  of  homogeneity  and  isotropy  [2]. 

3. 3  Derivation 

We  define  Dn(R)  as  the  refractive  index  spatial  covariance  function  where  we 
assume  that  the  random  refractive  index  field,  n(P),  is  statistically  homogenous  and 
isotropic  within  a  plane  that  is  parallel  to  the  tangential  plane  at  the  surface  of  a 
spherically-modeled  earth  for  the  optical  path  of  interest.  R  refers  to  the  Euclidean 
distance  between  two  points,  while  P  refers  to  a  point  in  three-space.  We  wish  to 
search  for  some  function,  Bn(R),  which  provides  the  best  statistical  description  of  a 
particular  set  of  Zernike  tilt  [24]  observations. 

3.3.1  Tilt  Covariance.  First,  we  define  the  covariance  operator  as 

cov[x, ?/]  =  E[xy\  -  E[x]E[y],  (3.1) 

and  we  define  unsealed  Zernike  tilt  as  the  quantity, 

a  -  JJ  <p[P±)P±  ■  xdPi,  (3.2) 

.4 

where  Pi  represents  a  point  in  two-space,  transverse  to  the  optical  (2)  axis,  and  ip 
represents  the  optical  phase  delay  between  two  points  in  space,  separated  by  distance 
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A.'  along  the  optical  axis  and  by  a  transverse  distance  of  \P±  (Az)  -  P± (0)|: 

Aj 

<P(P±)  A  !  n{P{z))dz  (3.3) 

o 

It  is  implied  by  using  rectilinearly  integrated  phase  delays  that  spatially  irradiance 
fluctuations  over  the  aperture  are  negligible  for  the  optical  system  under  consideration. 
For  a  more  detailed  analysis  of  this  assumption,  see  Appendix  A.  This  does  not 
preclude  its  application  to  conditions  of  strong  turbulence  as  defined  under  the  Rytov 
variance  criterion,  if  the  aperture  is  sufficiently  large.  Let  A  represent  the  aperture, 
and  let  the  plane  of  the  aperture  plane  be  defined  as  2  =  0,  while  the  object  plane 
is  defined  as  2  =  A2.  We  define  P(z)  to  be  a  3D  position  along  some  line  P  as  a 
function  of  position  along  the  optical  axis: 

m  =  {±(p„  -  a.) + ft.,  *(ft,  -p„)+p,.,x)  ,  (3.4) 

and  therefore  P±  is  related  to  P  by 

Pi  =  (Pro,  Pyo)  =  (P( 0)  •  P( 0)  •  V)  ■  (3.5) 

Explicitly  providing  an  argument,  implies  the  following: 

Pi(z)  =  (P(z)-x,P(z)-y).  (3.6) 

It  is  inferred  in  the  context  of  this  chapter  that  angle  brackets,  (  and  ).  are  used  to 
enclose  coordinates  for  a  point,  rather  than  to  denote  a  statistical  operation. 

We  declare  Bn(R)  to  be  the  refractive  index  covariance  function  under  homoge¬ 
nous  and  isotropic  turbulence  conditions.  Let  P\(z\)  represent  a  line  passing  through 
some  point  in  the  aperture  to  another  point  on  the  object  plane,  and  let  P2(z2)  repre¬ 
sent  a  second  line  passing  through  some  other  point  on  the  aperture  to  a  second  point 
of  interest  on  the  object  plane.  See  Figure  3.1.  Restrict  the  choice  of  Px  and  P>.  so 


41 


Figure  3.1:  Geometry  for  Tilt  Covariance. 
This  figure  depicts  the  geometry  for  computing 
tilt  covariance  for  two  separate  point  sources,  mea¬ 
sured  at  a  common  aperture. 


that  Pl(Az)  —  (— 4f,0, Ai)  while  P2{Az)  =  (-^,0, Ac).  Also,  define  the  following: 


Asi  =  |Pi(A,z)  —  Pi  (0)|, 

(3.7) 

As2  ^  \P2(Az)  -  P2(0)|. 

(3.8) 

Using  these  definitions,  we  can  express  the  tilt  covariance  between  two  points  on  the 
object  plane  separated  by  Ax  along  the  x-axis: 


Az 


Ba(Ax)  =k2  cov  ||  //ft  J  n{Px(z{))dzx  x  (Pu  -  x)dPu  J  , 

V.  A  0  , 

<  Az  ' 

Jf&J  n(p^))dz2  x  (PX2  •  x)d.P±2 

\  A  0  , 

2  As  A; 

=  ~Ez\  jj //^x  / y^cov[n(-Pi(2i)),n(P2(22))]  X  (3.9) 


a  a 


0  0 


X  (P±1  .  i)(Pi2 .  x)dPud.Pu 

2  Az  Az 

f:]  //  //AsiAs2X  /  /  Bn(\P>  -  P\\)dzxdz2 

A  A  0  0 

X  (Pi,  •  *)(Pi.2  •  x)dPudPl2, 


where  we  have  implied  the  dependency  of  on  the  aperture  function,  A\  the  optical 
path  length,  A z;  the  free-space  propagation  constant,  k;  and  any  parameters  which 
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might  pertain  to  the  refractive  index  covariance  function  Bn.  Let  us  also  imply  that 
for  a  sufficiently  long  optical  path,  we  may  translate  P\(Az)  and  P2(Az)  by  any 
amount  in  the  indirection  without  affecting  Ba( Ax),  so  long  as  the  points  are  still 
separated  by  Ax.  That  is  to  say,  we  assume  that  tilt  is  wide-sense  stationary  in  x. 
Approximating  (3.9)  for  any  sizable  path  by  using  six  nested  Riemann  sums  quickly 
becomes  intractable.  The  order  of  difficulty  grows  proportionately  to  the  square  of 
the  propagation  length,  Az.  In  the  next  subsection  we  will  analyze  a  more  practical 
numerical  approach  to  the  problem. 


3.3.2  Step  Interpolation.  It  turns  out  that  a  more  palatable  solution  results 
from  substituting  numerical  integration  with  a  combination  of  piecewise  symbolic 
integration  and  interpolation.  We  begin  by  defining  a  sequence  of  step  interpolates, 
{/*},  such  that 

Si  =  Bn(iAR),  (3.10) 

where  0  <  i  <  N  —  2.  Then,  we  approximate  Bn  for  positive  real  R,  using 

N-  2 

Bn(R)  «  J2  fi rect  (B  -  (2*  +  !))  ■  (3-11) 

i=0 

such  that  0  <  R.  <  (N  -  1)A R,  and  where 


rcct.(a:)  A  ) 


iAR  <  x  <  (i  +  1)A R 
elsewhere 


(3.12) 


After  substituting  the  step-interpolated  expression  for  B„  from  (3.11)  into  the  expres¬ 
sion  for  Be  from  (3.9),  the  innermost  double  integral  has  a  symbolic  solution  which 
depends  on  i,  the  index  of  the  corresponding  interpolate.  We  define  the  integral  /,  as 

Az  Az 

IdPi^Pi,,  Ax)  A  I  f  rect  _  (2-  +  1}\  (3.13) 

0  0 
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so  that  the  innermost  double  integral  from  (3.9)  becomes 


Az  A2 

I  -  j  I  B„{R  (zuz2,  Pu,P±2,Ax))  dzldz2 

0  0 
N- 2 

=  '£fiIi(P±l,Pu,  Ax). 

i=0 


(3-14) 


Since  P±,.  P±2,  and  Ax,  may  be  held  constant  with  respect  to  the  integration  defined 
in  (3.13)  ,  we  discard  them  from  our  notation  during  this  part  of  the  derivation.  Earlier, 
we  implicitly  defined  R{z\,z2)  to  be  the  Euclidean  distance  between  points  P\{z\)  and 
P2(z2).  This  quantity  is  formally  defined  here  as 


R{zi,z2)  -  [(21  ~z2)2 

+  [' mxizi  +  bXl  -  (mX2z2  +  bX2)f  (3.15) 

I 

+  [myi  z\  +  b;n  -  ( my2z2  +  bV2)f] 2 


where  mXj,  m,yi ,  ml2,  and  my2  are  the  respective  x-z  and  y-z  slopes  for  Pi  and 
P>  Similarly,  bXl,  byi ,  bX2,  and  bV2  are  the  x  and  y  intercepts  at  z  =  0.  Since 
R  is  a  real  scalar  function,  continuously  defined  over  the  real  z\-z2  plane,  there  is 
exactly  one  value  from  the  range  of  R.  corresponding  to  any  particular  point  from  the 
domain  of  R.  Therefore,  the  locus  of  points  satisfying  R  <  iAR  is  a  subset  of  the 
locus  of  points  satisfying  R  <  (i  +  1)A R.  We  define  5,  as  the  locus  of  points  from 
(21,22)  €  [0,  Ac]  x  [0.  As]  satisfying  R  <  iAR.  Then,  the  area  of  this  region,  A,,  is 
defined  as 


(3.10) 


Consequently,  we  can  form  the  equivalency, 


h  —  —  Aj. 


(3.17) 
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Next,  we  will  solve  for  A,  analytically.  Define  the  following  quantities  for  sub¬ 
stitution: 


°i  =  1  +ni2xi+m2yi: 

fl2  —  1  +  m2Xi  +  niy2 , 

an  =  2(1  +  mXlml2  +  myimy2), 

a  (3.18) 

ci  =  2  \mxi  (bXl  -bX2)  +  myi  (byi  -  by2)} , 

C2  =  2  [inX2  (bx ,  —  bX2)  +  mm  (byj  —  6y2)] , 

d  —  {bXl  —  bX2)2  -f  (by,  -  by.,)2  . 

Then,  rewrite  /?  as 


/?(2l,  22)  =[ai2j  +  Z2z\  -  «122l22  +  Ci2i  -  C222  +  rf] 2  .  (3.19) 


In  order  to  compute  A,,  we  must  first  determine  the  set  of  all  possible  points,  (zi,z2), 
such  that 

0  <  zi  <  A 2  (3.20) 

and 

0  <  22  <  A 2  (3.21) 


and 


Ol-sjf  +  z2z\  -  0122122  +  Cj2i  -  c222  +  d  <  V, 


(3.22) 


where 

u  =  (iAR)2. 


(3.23) 


(3.22)  comes  directly  from  the  inequality,  R  <  iAR.  Since  ci\  >  0,  we  can  rewrite  the 
expression  from  (3.22)  as  a  second-order  polynomial  in  21  whose  coefficients  depend 


Z\  +  2,  +  '^-c^+d-v  <  0 


(3.24) 
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To  complete  the  square,  we  define  b  and  h  such  that 


—  h  +  (21  +  b)1 


=  ^  +  (azSafl)*i  + 


a-i-j-c-yz-i+tt-v 

“1 


(3.25) 


This  leads  to  the  definitions 

rf-tA 
“>  )  ' 

Since  z\,  b ,  and  h  are  real,  (3.22)  becomes  simply 

-6  —  /ii  <  z\  <  -b  +  /is. 


I,  Aci-«X3*2 

‘  2a  i 


(3.2G) 


(3.27) 


We  introduce  another  set  of  substitutions  that  will  allow  us  to  cleanly  bound  z\  in 
terms  of  z2: 


% 


A 


ft  12 

'2a  j  > 


A  c? 

94  =  4^  “ 


ai 


(3.28) 


The  resulting  upper  and  lower  bounds  from  (3.27)  are 

-b  -  h-2  =q0z2  -qi-  [q2z\  -  q3z2  +  fc]*, 
-b  +  Jr-  =q0z2  -  qx  +  [q2z2  -  q3z2  +  q4]\ 


(3.29) 


which  implies  that 

q2zl  -  q3z2  +  q,\  >  0.  (3.30) 
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We  label 


P{z2)  =  q2zl  -  (hz2  +  q,i 


(3.31) 


as  the  characteristic  polynomial  for  z2.  From  (3.18)  and  (3.28)  it  is  clear  that  \q2\  <  1. 
Therefore,  the  special  case  where  q2  =  0  must  be  tested  prior  to  treating  (3.31)  as 
quadratic.  Next,  we  will  show  that  q2  <  0.  Expanding  q2  in  terms  of  its  defining 
slopes,  yields 


92  -i  (fe  -  “2) 

[(l+,,s:;;::rv2)  -  a + < + o 

(  m:r ,  )  ~f  (,ftj,  (  ^l/j  ~  Irl-y2  ) 

<0. 


(3.32) 


This  result  is  important  because  it  allows  us  to  directly  use  the  discriminant  corre¬ 
sponding  to  the  characteristic  polynomial  of  z2  to  confine  our  analysis  to  cases  where 
the  following  inequality  holds  true: 


qi  >  4 q2q-i- 


(3.33) 


All  other  cases  are  not  physically  realizable. 

The  next  step  in  computing  A,;  is  to  identify  the  allowable  interval  for  z2.  Since 
q2  <  0,  the  characteristic  polynomial  for  z2  is  concave  down.  For  cases  where  (3.33) 
is  true,  we  can  say  that  there  is  exactly  one  or  zero  such  intervals  on  z2  contributing 
to  A,  satisfying  both  (3.30)  and  (3.21).  We  label  the  upper  and  lower  bounds  of  this 
interval,  C2  and  C->+>  respectively.  If  for  the  particular  case  of  interest  there  are  no 
such  intervals,  then  A*  =  0.  We  continue  our  analysis,  assuming  that  the  interval 
[C~  -  C+]  exists  for  the  given  choice  of  i. 

It  may  be  necessary  to  partition  [C2",  C2I  into  multiple  mutually  exclusive  subin¬ 
tervals  for  which  the  limits  on  Z\  satisfying  (3.20)  arid  (3.22)  are  fixed.  In  other  words, 
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the  limits  of  integration  for  z\  may  change,  depending  on  the  value  of  z2.  We  rewrite 
(3.16)  to  reflect  the  use  of  such  partitions: 


A, 


A  —0i  J  I  dz\dz2  +  02  I  j 

Si  o  s2  o 


-b+Vh  A  z  -b+Vh 

dz\dz2  +  0:s  J  j  dzldz2  +  0A  J  j  dzxdz2 


^3  -b—t/il 


S.i 


— /3i  AS] As  + /?2  j b  +  h- )dz2  +  fyAz  +  03  j (b  +  h^)dz2  +  20A  J  h*dz2, 
s2  s3  s4 


(3.34) 


where  AE,  =  fs^  dz2.  A  boolean  (0)  coefficient  is  placed  before  each  of  the  integrals 
corresponding  to  the  four  possible  partitions.  Each  0  coefficient  is  restricted  to  have 
a  value  of  either  i  or  0.  We  will  reduce  (3.34)  by  defining  the  definite  integrals: 


and 


(3.35) 


c2+ 

Q/.(C2-,C 2+)=  J  h? dz2 

C2 


cl 

J  [9222  ~  Q3Z2  +  <74] 5  dz2. 
^2 


(3.36) 
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If  we  presume  that  q2  =  0,  then  (3.36)  evaluates  to  [33] 


Qh  = 


2(?4  +  93 


\  3 


'<2+ 


3?3 


(3.37) 


Otherwise,  for  q2  <  0,  we  have  [33] 


Qh 


-  Kz*  +  71]2  (f  -  i|) 

1m  1 2r/!  s2  -  <m~J  +  2[g2z£  -  qzz2 
»</■?  L 


(3.38) 


Since  the  integrand  is  known  to  be  real  over  the  entire  interval  of  integration,  we 
can  safely  analyze  just  the  real  part.  This  is  computationally  more  efficient  because 
Qh(0,  Q )  and  Qh[ 0,C^)  are  complex  conjugates.  The  resulting  expression  is 


Qh  = 


'  '  Sn" 


1 Q 


x  axctan 


f  2|92|-^2  +  ^3 1 <72 1  - 
(2[q2z%  -  q3z2  +  /L- 


Using  these  integral  function  definitions,  we  rewrite  (3.34)  as 


Ai  =AAE,A  z  +  ft[-Q*(£a)  +  Qh(  E2)] 

+  #j[Az  +  C?ft(S.3)  +  Qft(S3)] 

+  2/34Q/l(£4). 


(3.39) 


(3.40) 


The  mechanical  details  on  how  to  solve  for  the  exact  boundaries  of  the  z2  partitions, 
S|  -S4,  are  beyond  the  scope  of  this  document.  Suffice  it  to  say  that  eacli  of  the  par- 
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tition  boundaries  has  an  analytical  solution  whose  resulting  form  may  he  determined 
by  evaluating  a  set  of  test  conditions  in  a  particular  order. 

At  this  point,  we  rewrite  (3.9)  as 

N- 2 

As,  As2  J2  /«(A+ l  -  A)  (Pu  •  *)  (P±2  ■  ®)  dP_L,dPu, 
1=0 

(3.41) 

and  observe  that  the  two  innermost  integrals  from  (3.9)  have  been  condensed  into  a 
single  sum  rather  than  into  two  nested  Riemann  sums.  Furthermore,  the  number  of 
elements  in  tliis  sum  can  remain  fixed  with  respect  to  the  optical  path  length,  A z, 
without  further  loss  of  fidelity.  Instead,  the  complexity  of  this  summation  depends 
on  the  values  of  l0  and  L0.  To  see  why  this  is  true,  we  observe  from  (3.11)  that  for 
R.  >  (N  —  1)A R,  Bn  is  approximated  as  zero.  ( N  —  1)A R  should  be  set  to  a  value 
greater  than  or  equal  to  L;.  Additionally,  the  resolution  of  Bn  is  limited  to  the  sample 
spacing,  A R.  The  exact  value  used  for  A/?  should  be  chosen  based  on  /0.  Neither  of 
these  quantities  depends  on  path  length  or  geometry.  Once  A R  and  N  are  fixed,  it  is 
possible  to  tabulate  Bn( Ax)  for  any  given  path  and  imaging  system  and  for  a  desired 
set  of  point  separations,  {Ax*}. 

3.3.3  Spectral  Decomposition  of  Refractive  Index  Covariance  Function- 
While  the  previous  expression  for  tilt  covariance  in  (3.41)  may  have  some  utility  in 
determining  the  characteristic  parameters  of  a  particular  refractive  index  covariance 
function,  the  results  may  be  limited  to  a  particular  set  of  atmospheric  conditions, 
and  to  a  particular  path  of  observation.  In  this  section,  we  will  derive  an  efficient 
means  of  spectral  tabulation  that  facilitates  fast  computation  of  tilt  covariance  for 
infinitely  many  refractive  index  covariance  functions.  The  traditional  starting  point 
for  statistical  analysis  of  refractive  index  fluctuations  is  in  the  spectral  domain.  In 
this  derivation,  we  conversely  start  in  the  spatial  domain,  since  our  observable  data 
(tilt  covariance)  are  most  closely  related  to  that  domain. 


Ba{  Ax)  = 


MU  II 
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A  complex-valued  discrete  signal,  x,  is  related  to  its  discrete  spectrum,  X,  by 
the  following  pair  of  relationships:  [19] 

N 

Xk2  =  **,  exp  { -j g(fci  -  l)(/c2  -  1) }  (3.42) 

*1=1 

and 

1  N " 

Xb  exp  {|$S  (*,  -  Ij (fe  -  1)}  .  (3.43) 

S  *2=1 

In  order  to  take  advantage  of  the  discrete  Fourier  transform  (DFT)  relationship,  we 
must  somehow  sample  our  continuous  representation  of  Dn.  One  trivial  approach  is 
to  apply  a  uniform  sampling  grid  with  samples  spaced  at  A R.  However,  this  does  not 
bode  well  with  cases  where  lo  <C  Lq,  and  such  cases  are  expected  to  be  the  norm.  We 
explore  an  alternative,  where  we  can  reduce  the  total  number  of  samples,  N„.  This  is 
accomplished  by  transforming  the  /Taxis  prior  to  applying  a  uniform  sampling  grid 
and  executing  the  DFT.  For  that  purpose,  we  define  a  function,  </(•),  and  its  inverse, 
</-i  (•),  such  that  our  original  spatial  variable,  R,  is  transformed  into  Rg: 

Ra  = 9{R ), 

R=9~\Rg), 

We  require  that  </(■)  be  monotonically  increasing.  Within  the  ^-transformed  space, 
we  define  the  uniform  samphng  distance  to  be  A Rg.  For  this  particular  problem, 
we  select  g(-)  using  the  following  criteria: 


and 


5(0)  =  0, 


5  (£o„„,  J  —  L( i 


(3.45) 


Bnia-'ihARg))  -  Bnig-'dk!  -  l)ARg))\  «  K,  (3.46) 
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where  L0max  is  the  maximum  outer  scale  that  can  be  modeled  by  the  resulting  spectral 
decomposition.  The  third  criterion  ensures  that  each  transformed  sample  interval 
represents  approximately  equal  change  in  the  value  of  Bn(R).  This  criterion  allows  us 
to  efficiently  sample  a  refractive  index  covariance  function  with  near-zero  inner  scale 
without  significant  aliasing.  In  order  to  solve  for  g(-)  we  need  to  assume  some  initial 
form  for  Bn.  For  this  work,  we  have  chosen  to  use  Bn  resulting  from  the  von  Karnian 
spectrum  with  l.0  =  0.  It  can  be  empirically  demonstrated  that  the  von  Karman 
refractive  index  covariance  function  can  be  closely  approximated  by 

Bn.,0(fl)»CM“677/n(g),  (3.47) 

for  0  <  R  <  oo  and  /0  =  0,  where 


f  (D\  —  _ 2.SH.'15 _ 

(/;0.7023  +2.1,195R1.0.157  +  1  2355)5 


(3.48) 


Once  g( ■ )  has  been  determined,  we  can  construct  the  transformed  spectral  basis  set. 
We  will  use  linear  interpolation  to  define  the  spectral  basis  continuously  in  R ,  in 
relation  to  the  discrete  representation  in  the  transformed  space,  Rg.  For  a  particular 
value  of  R. :  0  <  R  <  we  define  the  nearest  neighboring  samples  corresponding 
to  the  transformed  space  as  R„  and  Ri,: 


R n  =g~l 
Rb  — 3  1 


A  Rv 

w 

A  Ra 


(3.49) 
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Relating  (3.49)  and  the  complex  exponentials  of  (3.43),  we  define  the  following  angular 
quantities  of  substitution: 


"«  -  1  )(k2  ~  1) 

-  1), 

**±£(*l)(fe-l) 

- 1). 


(3.50) 


where  we  have  implied  (Aq  -  l)Ai?g  =  Rn  and  A/?,,  =  /4,  such  that  k\  is  dependant 
on  /?  by  (3.49).  These  assignments  provide  for  the  following  continuous  definitions  for 
linearly  interpolated  real  and  imaginary  spectral  basis  functions  used  to  decompose 
Bn(R): 


(R)  [fcfe  (C0SK)  ~  cos(w6))  +  cos(cua)]  , 

W^k2(R)  =  -  i  [it^;  (™k)  ~  sin(w6))  +  sin(wa)]  . 

The  bases  in  (3.51)  come  directly  from  applying  (3.50)  to  (3.43).  W3hi(R)  is  the 
real  component  of  the  kfj'  spectral  basis  function,  and  M/o*2(H)  is  the  corresponding 
imaginary  component.  According  to  these  definitions,  the  continuous  function  Bn(R) 
may  be  approximately  reconstructed  using  the  spectral  decomposition, 


N 


Bn(R)  =  E  (*  {®»*  }  W*. (R)  +  9  {*„„  }  W*. (fl>) 

fc2=l 

*Bn(R), 


(3.52) 


where  the  kf  complex  spectral  component  is  defined  from  (3.42)  as 

AT, 

^  *£*  tr'tfc  -  »%))  exp  {-if  (*.  -  life  -  1)}  -  (3.53) 

ki=\ 

By  applying  (3.41)  to  each  of  the  Ns  real  and  imaginary  spectral  basis  functions 
defined  in  (3.51),  a  table  can  be  generated  for  Bn( Ax)  which  can  be  used  on  any  form 
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of  B„(R),  in  the  following  manner: 


Ba(Ax)  =  XI  {SatoM Ax)3?  {$„fc2 }  +  BoMz(Ax)Z  {<f>nfc2 })  .  (3.54) 

fcs= 1 

Ba,k2,n{ Ax)  and  5„i*;.2icJ(A.t)  are  the  tabulated  tilt  covariance  values  associated  with 
a  source  separation  of  Arc  and  the  k^1  real  and  imaginary  spectral  bases. 


3.3-4  3-Parameter  von  Karman  Model  for  Bn.  Next,  we  introduce  an  effi¬ 
cient  method  for  tabulating  the  refractive  index  covariance  function  associated  with 
the  von  Karman  spectrum.  The  von  Karman  spectrum  for  refractive  index  fluctua¬ 
tions  is  given  by  [2] 


<hn,»  =  0.033C2 


K2  + 


(*)' 


exp  {— 0.02853k.2 /o} 


(3.55) 


where  the  k  denotes  one-dimensional  spatial  frequency.  Bn(R )  and  $„(«;)  are  related 

by  [2] 

oc 

47T  f 

Bn(R)  =  —  /  <!>„(«)  sin(/c/?)KdK.  (3.56) 

0 

It  can  be  empirically  shown  that  the  von  Karman  model  very  nearly  obeys  the  fol¬ 
lowing  compound  scaling  property: 


Bnvk(C2n,l0,L0,R) 


•‘1  r  0.6677 


CznL0 


n  ( 1  2il  i  JL\ 
^  Lo'  ’  Lo) 


(3.57) 


This  property  allows  us  to  tabulate  Bn  as  a  function  of  only  R  and  /n.  For  our 
purposes,  we  used  a  1000-point  log-spaced  sampling  grid  in  R  and  a  101-point  uniform 
grid  for  A  g  [().  0.01  j.  Linear  interpolation  is  used  to  obtain  Bn  for  other  combinations 
of  R  and  Iq.  This  approach  is  fast  and  convenient.  (3.57)  is  a  valid  approximation  for 
values  of  L0  on  the  interval  [1, 100].  For  this  work,  units  of  meters  were  applied  to  all 
spatial  quantities. 
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3-4  Experiment 

We  designed  a  simple  experiment  to  demonstrate  liow  a  simple  low-cost  ap¬ 
paratus  can  be  deployed  to  make  field  measurements  of  the  refractive  index  spatial 
covariance  function,  Bn. 

3J,.l  Description.  The  imaging  system  consisted  of  a  Photometrics®512B 
high-speed  broadband  CCD  camera  [26],  coupled  to  a  Meade®8-in  LX200GPS  Scmidt- 
Cassegrain  telescope  [20] .  The  LX200GPS  Inis  a  2  m  focal  length,  and  was  custom- 
fitted  with  a  4  cm  circular  aperture  stop.  Camera  specifications  include  a  512  x  512 
detector  array  of  16  pm  x  16  pm  broadband  pixels.  Peak  pixel  quantum  efficiency  is 
achieved  at  a  wavelength  of  approximately  575  nm,  with  quantum  efficiency  ranging 
0.90  <  i]  <  0.93,  for  500  nm  <  A  <  675  nm.  A  532  nm  narrowband  filter  was  placed 
directly  over  the  aperture  stop  to  reduce  the  photon  flux  at  the  receiver.  The  camera 
transmits  uncompressed  image  data  to  a  PC  via  a  PCI  card  at  a  rate  of  28  fps.  A 
photograph  of  this  apparatus  is  shown  in  Figure  3.2. 

The  imaging  apparatus  was  configured  to  view  a  passive  array  of  extended 
sources  through  a  100  m  atmospheric  path  at  1.25  m  above  ground.  The  source  array 
was  configured  to  be  one  dimensional,  parallel  to  the  ground,  and  perpendicular  to 
the  optical  axis.  Figure  3.3  shows  the  pattern  for  the  source  array.  The  distribution 
of  sources  was  stochastically  optimized  to  achieve  the  maximum  number  of  distinct 
separation  distances  between  source  pairs  that  are  uniformly-distributed  over  22.5  cm. 
Table  3.1  gives  the  relative  positions  of  each  source,  while  Figure  3.4  shows  thirty-three 
corresponding  usable  distances  that  follow  a  nearly-uniform  distribution.  The  reason 
for  favoring  a  uniform  distribution  of  source-pair  distances  is  that  the  distribution  of 
relevant  information  with  respect  to  Ax  is  unknown.  Therefore,  this  is  the  distribution 
which  will  most  likely  maximize  the  total  information  content  in  the  data  collected. 
Likewise,  the  reason  for  maximizing  the  number  of  distinct  source-pair  distances  is 
also  to  maximize  the  total  amount  of  useful  information  recorded.  The  span  of  22.5  cm 
is  designed  to  nearly  fill  the  field  of  view  for  an  optical  path  length  (OPL)  of  100  m, 
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Figure  3.2:  Imaging  Apparatus.  This  figure  shows  the  setup  for  the  imaging  apparatus 
used  in  the  experiment. 

while  still  allowing  some  room  for  relative  motion  between  the  sources  and  the  imaging 
instrument. 


3.4.2  Observables.  In  order  to  account  for  such  relative  platform  motion, 
the  observed  tilt  covariance,  D&,  is  defined  as 


Bc>  =  cov[d, ,  d'i], 


(3.58) 


where, 


di  =  (*1  +  av 
a2  =  a2  +  a,,. 


(3.59) 
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Figure  3.3:  Passive  Source  Array.  This  figure  shows  pattern  used  for  the  source  array. 
There  are  ten  sources  placed  in  a  manner  which  gives  thirty-three  uniformly-distributed 
distances  of  separation. 


Table  3.1:  Relative  Source  Positions 


Source 

Position  (cm) 

1 

0.0000 

2 

1.5000 

3 

3.0000 

4 

4.7539 

5 

8.6730 

6 

13.0252 

7 

18.0000 

8 

19.5000 

9 

20.1000 

10 

22.5000 

Consequently, 

Ba  =  Ba  +  vax[a„],  (3.60) 

where  a,t  is  the  noise  in  tilt  due  to  translational  platform  motion.  The  platform 
motion  is  assumed  to  be  statistically  independent  from  atmospheric  tilt.  Pitch  and 
roll  do  not  affect  the  observations,  and  yaw  is  assumed  to  be  negligible.  var[cv,;]  is  a 
fourth  unwanted  parameter  that  vail  be  estimated  along  with  the  other  turbulence 
parameters.  The  random  platform  motion,  av,  is  constant  with  respect  to  source 
position  for  simultaneous  observations:  =  av(xi )  =  ar,(.r2). 

Next,  we  address  spatial  aliasing  of  the  system.  Using  Nyquist  theory,  the 
aliasing  criterion  for  this  system  is  AJllx  <  We  get  this  relationship  from  the  fact 
that  the  maximum  spatial  frequency  of  the  electromagnetic  (EM)  field  at  the  detector 
is  ^  111].  The  intensity  pattern  has  twice  the  maximum  spatial  frequency  as  the  EM 
field,  and  Nyquist  sampling  requires  the  maximum  sampling  distance  to  be  twice  the 
maximum  frequency  of  the  intensity  pattern.  For  our  particular  imaging  system,  we 
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Nearly  Uniform  Distribution  of  Source  Separations 


Figure  3.4:  Distribution  of  Source  Separations.  This  plot  displays  the  distribution  of  distances 
between  pairs  of  sources,  from  Figure  3.3.  The  distances  are  stochastically  optimized  to  be  maximal 
in  number  and  nearly  uniform  in  distribution. 

would  need  a  pixel  pitch  of  A  pi*  <  6.65  pm  to  eliminate  all  aliasing.  Our  design  is 
therefore  subject  to  some  aliasing.  However,  the  spatial  frequency  power  spectrum 
tapers  off  linearly,  since  it  is  sealed  by  the  autocorrelation  of  the  pupil  function. 
Using  this  analysis,  we  can  directly  compute  the  maximum  total  power  of  aliasing,  as 
a  fraction  of  the  maximum  total  spatial  power  in  the  image.  The  maximum  spatial 
image  power  is  proportionate  to 


K, 
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CC 


3A2  ' 

nyq 


(3.01) 
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where  Anw  is  the  maximum  pixel  pitch  allowed  by  the  Nyquist  criterion.  This  defini¬ 
tion  comes  directly  from  using  the  formula  for  the  volume  of  a  right-circular  cone  [28]. 
The  amount  of  non-aliased  spatial  image  power  for  the  image  system  we  used  is  pro¬ 
portionate  to 


3A^ 


1  _  Any(/ 
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pix 
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A2-  ’ 
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(3.62) 


which  comes  from  computing  the  volume  of  the  same  cone,  intersected  with  a  right- 
circular  cylinder.  The  resulting  ratio  of  aliased  spatial  power  to  total  spatial  image 
power  is  given  by 


5=1- 
=1  - 


(3.63) 


For  our  particular  instrument,  the  fraction  of  aliased  power  for  the  worst-ease  image 
scenario  is  5  =  0.6254.  Because  this  is  quite  poor,  we  selected  rather  large  (1  cm) 
diameter  sources.  The  extended  sources  enable  us  to  more  accurately  estimate  tilt 
because  they  produce  images  with  lower  overall  spatial  frequency  content,  while  si¬ 
multaneously  improving  the  signal-to-noise  ratio  (SNR)  of  the  system.  We  purposely 
did  not  further  reduce  the  aperture  of  the  system  to  avoid  all  aliasing.  Doing  so  would 
have  degraded  the  resolving  power  of  the  imaging  system.  This  would  have  required 
larger,  and  hence  fewer  sources,  due  to  optical  turbulence  effects.  As  the  aperture 
becomes  smaller,  the  potential  for  scintillation  to  occur  is  greater.  Using  a  larger, 
slightly  aliased  aperture  allows  us  to  use  less  robust  image  processing  techniques  that 
do  not  need  to  account,  for  scintillation.  Alternatively,  we  could  have  selected  a  differ¬ 
ent  camera  with  smaller  pixels.  However,  this  would  require  lengthening  the  exposure 
time,  in  order  to  achieve  sufficient  photon  flux.  This  is  may  not  be  desirable  because 
the  atmosphere  is  continually  changing.  Short  exposure  times  allow  us  to  assume  the 
atmosphere  is  frozen  during  the  integration  period  for  each  frame.  While  our  system 
is  not  parametrically  optimized,  it  is  a  sufficient  prototype  for  making  measurements 
of  tilt. 
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The  image  processing  algorithm  to  extract  tilt  observations  from  a  sequence  of 
images  is  straightforward.  First,  the  horizontal  and  vertical  profiles  are  computed  for 
each  frame,  by  summing  all  of  the  columns  and  then  all  of  the  rows.  Then,  peaks  are 
detected  from  the  horizontal  profile,  detecting  the  horizontal  location  of  each  source 
image.  Next,  the  vertical  location  of  the  row  of  sources  is  detected  from  the  horizontal 
profile.  At  this  point,  a  small  square  mask  is  used  to  extract  an  isolated  image  of  each 
source.  The  centroid  is  computed  for  each  source,  as  a  preliminary  estimate  for  the 
source’s  location.  Then,  each  source  is  extracted  once  more  from  the  original  image, 
using  a  new  square,  centered  at  the  centroid.  Next,  a  horizontal  profile  is  computed 
for  each  source  image.  Each  profile  is  interpolated  by  a  factor  of  ten,  before  it  is 
correlated  with  a  truncated  normal  probability  density  function  (pdf).  The  peak 
correlation  is  used  as  the  subpixel  estimate  for  the  location  of  the  image  source.  Now, 
we  will  show  that  the  resulting  quantity  is  proportionate  to  tilt  if  momentarily  neglect 
the  noise  component. 

Consider  the  diagram  in  Figure  3.5,  which  relates  angle  of  arrival  to  diffraction- 
less  image  displacement,  as  predicted  by  free-space  geometric  optics.  We  will  assume 
that  the  refractive  index  is  constant  with  respect  to  wavelength,  for  all  media  under 
consideration.  Geometric  tilt,  ag,  is  closely  related  to  the  angle  of  arrival,  9.  For  our 
application,  we  have  the  condition,  |Ax;|  -C  |/|,  which  leads  to  the  approximation, 
sin(i9)  =  All  ~  0.  We  formally  define  geometric  tilt  as  the  optical  phase  delay  per 
unit  distance  across  the  aperture  in  the  x-direction: 
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Figure  3.5:  Geometric  Tilt.  This  diagram  depicts  the  relationship  between  angle  of  arrival  and 
image  displacement  from  the  optical  axis. 

where  k  =  ^  is  the  propagation  constant  for  free  space.  It  follows  from  our  definition 
of  geometric  tilt,  that  the  phase  over  the  aperture  due  to  ag  is  given  by 


ip{x,  y)  =  xag.  (3.65) 


From  (3.2),  the  corresponding  amount  of  unsealed  Zernike  tilt  is 


(3.66) 
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Transforming  from  rectangular  to  polar  coordinates,  by  making  the  substitutions, 
dxdy  =  rdrdd  and  x  =  rcos(0),  yields 
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where  ?e  is  a.  rational  number  representing  the  number  of  pixels  by  which  the  image  is 
displaced  in  the  ^-direction.  This  shows  that  instantaneous  tilt  observations  can  be 
obtained  by  scaling  corresponding  observations  of  instantaneous  image  displacement 
due  to  each  source  in  the  ^-direction. 


3.4-3  Tabulation  and  Estimation.  We  introduce  the  following  maximum- 
likelihood  estimator  for  covariance  of  observed  tilt,  Z/y,  assuming  that  tilt  is  normally 
distributed.  Given  m  independent  observations  of  tilt  for  two  sources,  located  at  X\ 
and  x-2,  such  that  Ax  =  |.ti  —  x-z\,  we  define  the  estimator  for  tilt  observed  covariance 
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where  cvcii ^  is  the  observed  tilt  due  to  source  1  from  the  Ith  image  frame,  and  <*2,/  is  the 
corresponding  observation  of  tilt  due  to  source  2.  For  a  particular  set  of  estimated 
turbulence  parameters,  we  define  the  corresponding  tilt  covariance  to  be  Ba(Ax), 
which  is  computed  using  the  pretabulated  coefficients  given  by  (3.54).  We  employ  the 
implementation  of  the  Nelder-Mead  simplex  algorithm  packaged  with  MATLAB®  [19] 
to  search  for  the  values  of  /<,  and  Lo  which  most  closely  describe  the  observed  data.  For 
any  given  choice  of  /0  and  L0,  we  solve  for  the  unconstrained  least  squares  estimates 
of  Cl  and  var[cvT,],  treating  them  as  gain  and  bias,  respectively.  The  cost  function 
used  for  this  process  is  normalized  mean  squared  error  in  zero-mean  tilt  covariance: 


771—1  ^ 

E  BUAxi) 

i=0 

where  <5^  «  var[«,,]  is  the  estimated  motion  variance,  and  I30vk(A:rt)  is  the  tilt 
covariance  predicted  by  the  von  Karman  model. 

3./,.4  Experimental  Conditions.  The  experiment  was  performed  on  June  15, 
2006,  beginning  at  1544  EDT  and  lasted  for  20  minutes.  The  optical  path  measured 
100  m,  and  the  path  was  1.25  m  above  ground.  The  surface  was  asphalt,  subject  to 
direct  solar  radiation,  at  an  equilibrium  temperature  of  47.22  C.  The  skv  was  clear, 
providing  for  irradiance  conditions  that  were  approximately  temporally  invariant. 
The  prevailing  cross-wind  measured  2.0  m/s  with  gusts  up  to  2.5  m/s.  The  mean  air 
temperature  was  28.6  C.  Relative  humidity  was  at  27%,  and  the  atmospheric  pressure 
was  stable  at  1028  rnbars.  The  operating  wavelength  was  A  =  532  run. 

3. 5  Results 

Atmospheric  turbulence  is  a  nonstationary  process.  However,  we  treat  tilt  as 
locally  stationary  over  short  segments  of  time.  Figure  3.6  shows  a  histogram  of  differ- 
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ential  tilt  observations  extracted  from  an  arbitrary  block  of  240  s.  Another  important 


Comparison  of  Differential  Tilt  Distribution  to  Normal  Distribution 


<x(0)  -  a(22.5  cm) 


Figure  3.6:  Differential  Tilt  Distribution.  This  figure  confirms  that  the  subset  of  data  chosen 
tor  analysis  is  approximately  normally  distributed. 

point  to  consider  is  that  the  data  is  temporally  correlated.  Ergodicity  would  imply 
that  discrete-time  statistics  are  equivalent  to  ensemble  statistics.  For  a  discrete¬ 
time  signal  that  is  ergodic,  we  could  approximate  ensemble  statistics  by  computing 
temporal  statistics  over  a  sufficiently  long  interval.  Unfortunately,  our  signal  is  not 
ergodic,  because  it  is  nonstationary.  Instead,  we  loosely  assume  that  first  and  second- 
order  ensemble  statistics  may  be  approximated  by  estimating  the  corresponding  time 
statistics  from  a  sufficiently  long  interval.  Clearly,  any  data  segment  may  be  randomly 
reordered  without  affecting  its  ensemble  statistics.  The  new  permutation  would  result 
in  a  temporally  uncorrelated  sequence,  whose  members  can  be  treated  as  independent, 
observations.  If  the  observation  period  were  sufficiently  long,  then  a  histogram  of  the 
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data  would  appear  unimodal  and  symmetric.  Since  Figure  3.G  shows  exactly  that, 
we  have  no  evidence  that  our  assumption  is  faulty.  While,  we  have  not  validated  our 
initial  assumption,  lack  of  any  disproof  is  sufficient  reason  to  continue. 

Because  we  did  not  know  exactly  what  value  to  use  for  L0wax  in  generating 
our  table  of  spectral  tilt  coefficients,  we  produced  several  tables  corresponding  to 
Ln„,„.  €  {1.0  in,  2.5  m,  5  m,  10  m, 40  m,  100  m}.  Each  table  was  produced,  using  a  PC 
equipped  with  dual  2.4  GHz  Xeon©processors  with  hyper-threading  enabled.  All  four 
virtual  processors  cooperatively  produce  a  single  table  in  approximately  106  hours. 
Once  produced,  the  tables  may  be  reused  for  the  same  experiment  under  different 
turbulence  conditions.  One  strategy  for  selecting  the  correct  table  is  to  start  with 
the  table  corresponding  to  the  largest  L0mas  and  to  decrease  L0lliax  to  the  smallest 
tabulated  value  which  is  greater  than  the  estimated  value  of  L0. 

Figure  3.7  shows  the  results  from  estimating  Bn,  using  L0  =  40  m.  The  optimal 
solution  for  Bn  achieved  by  the  von  Karman  covariance  model  for  the  whole  data  set 
is  given  in  Figure  3.8.  The  parameters  corresponding  to  the  optimal  solution  for  the 
ensemble  of  data  are  as  follows: 

Cl  5=6.01  •  10_9m'l, 
lo  5=17.9  mm, 

L0  «15.5  m. 

By  contrast,  the  average  parameters  corresponding  to  the  ensemble  of  all  240  s  win¬ 
dows  shown  in  Figure  3.7  are 

Cl  5=6.10-  10-9m~t, 
l0  5=13.6  mm, 

L0  5x14.8  m. 

The  two  sets  of  results  appear  to  be  consistent.  Figure  3.9  compares  the  estimated 
tilt  covariance  for  the  entire  data  set  to  the  actual  observable  data.  The  data  includes 
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Temporal  Variations  in  Turbulence  Parameter  Estimates 
(t  „  =  240  s,  At  =  10  s) 

'  window 
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Figure  3.7:  Windowed  Estimates  of  Turbulence  Parameters.  This  figure  displays  the  windowed 
estimates  of  /o,  Lo,  and  C\  corresponding  to  B„(R),  optimized  for  Lomax  =  40  m. 

the  effects  of  platform  motion.  According  to  the  figure,  the  von  Karman  models  the 
observed  data  very  closely.  Much  prior  literature  has  been  dedicated  to  presenting 
theory  and  results  based  upon  the  much  simpler  Kolmogorov  model.  It  has  been 
unknown  until  now,  how  well  the  |-Power  Law  describes  refractive  index  fluctuations 
in  strong  turbulence  conditions  near  the  ground.  Figure  3.10  provides  soiue  initial 
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Optimal  Solution  (of  Bo{R) 


Figure  3.8:  Optimal  von  Karman  Covariance  Function.  This  figure  displays  the  von  Kilrman 
Bn  that  optimally  describes  the  observed  data. 

evidence  to  debunk  the  use  of  Kolmogorov  theory  for  terrestrial  imaging  applications, 
especially  where  optical  turbulence  is  likely  to  be  strong. 

3. 6  Conclusions 

This  research  could  significantly  impact  the  development  of  future  applications 
in  terrestrial  remote  sensing.  We  have  presented  a  low-cost  method  which  allows 
the  refractive  index  covariance  function  to  be  passively  measured ,  assuming  a  von 
Karman  model.  The  method  does  not  rely  on  scintillation  phenomenology  and  does 
not  require  the  traditional  paraxial  assumption  to  be  used.  While  the  current  theory 
depends  on  statistically  homogenous  and  isotropic  coplanar  conditions,  there  may 
be  immediate  applications  where  such  conditions  are  guaranteed.  Certainly,  this 
measurement  technique  offers  a  potential  improvement  over  existing  Kolmogorov- 
based  techniques  for  optical  turbulence  measurement  and  characterization.  Future 
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Comparison  of  Tilt  Covariance  Model  to  Observations 
Normalized  MSE:  3.651  le-005 


Figure  3.9:  Tilt  Covariance.  This  figure  compares  the  tilt  covariance  estimated  directly  from 
all  observations  to  the  tilt  covariance  that  corresponds  to  Bn  from  Figure  3.8. 

work  will  be  directed  toward  establishing  specific  criteria  for  which  scintillation  effects 
can  be  considered  negligible.  The  work  presented  in  this  chapter  was  co-sponsored  by 
the  Air  Force  Institute  of  Technology  (AFIT)  Department  of  Electrical  Engineering 
(ENG)  and  by  the  Air  Force  Research  Laboratory  (AFRL)  Electro-Optical  Combat 
Identification  Technology  Branch  (SNJM). 
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4.  Summary 


This  dissertation  documented  three  new  contributions  to  help  advance  the  state 
of  the  art  in  simulation,  modeling,  and  estimation  of  optical  turbulence.  The 
underlying  theme  in  this  research  was  to  revisit  historical  assumptions  applied  to  pre¬ 
viously  intractable  computations  related  to  this  class  of  problems.  In  Chapter  2,  we 
presented  the  two  contributions.  The  first  contribution  is  the  formalization,  optimiza¬ 
tion.  and  validation  of  a  modeling  technique  called  successively  conditioned  rendering 
(SCR).  The  second  contribution  is  the  derivation  of  the  covariance  delineation  the¬ 
orem,  which  provides  theoretical  bounds  on  the  error  associated  with  successively 
conditioned  rendering.  It  was  shown  how  the  covariance  delineation  theorem  can  be 
used  as  a  computationally  efficient  tool  for  optimizing  the  performance  of  SCR  al¬ 
gorithms.  In  our  particular  example,  we  compared  two  different  SCR  algorithms  for 
simulating  random  phase  screens.  The  theorem  correctly  predicted  which  algorithm 
would  perform  more  accurately. 

Chapter  3  documented  a  third  contribution,  by  introducing  a  new  method  for 
estimating  optical  turbulence  parameters,  and  demonstrating  the  method  using  exper¬ 
imental  data.  For  our  particular  experiment  we  passively  imaged  an  array  of  sources 
through  a  100  m  path  at  1.25  m  above  the  solar-heated  tarmac  at  28  fps.  We  used 
10,000  samples  of  data  assumed  to  be  locally  stationary  for  parameter  estimation.  As¬ 
suming  a  von  Karman  turbulence  model,  we  estimated  the  atmospheric  parameters 
for  these  conditions  to  be  C%  «  6.10  •  10-9  m~t,  l0  «  13.6  mm,  and  L0  ~  14.8  m. 

4-1  Value  to  the  Air  Force 

Collectively,  these  three  contributions  are  suitable  for  use  in  several  Air  Force 
applications.  Contributions  1  and  2  provide  a  methodology  for  optimizing  SCR  al¬ 
gorithms.  By  applying  SCR.  to  phase  screen  simulation,  Air  Force  systems  engineers 
can  more  efficiently  and  accurately  characterize  optical  sensor  performance.  A  few 
examples  are  laser  vibrornetry,  imaging  systems,  and  deconvolution  algorithms. 
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Contribution  3  is  a  turnkey  method  for  realtime  estimation  of  optical  turbulence 
parameters.  This  method  can  be  used  to  support  field  tests,  and  it  can  also  be  used  to 
characterize  fade  environments  for  operational  laser  communication  systems.  Iligh- 
energy  laser  weapons  may  use  this  method  for  estimating  turbulence  parameters  by 
projecting  an  array  of  low-energy  laser  illuminators  onto  a  target.  The  imagery  of  the 
spot  array  could  also  be  used  to  solve  for  the  range  and  orientation  of  the  target,  given 
the  target,  geometry.  For  scientific  monitoring,  daytime  C%  profiles  may  be  obtained 
by  placing  a  low-weight  variant  of  the  instrument  on  a  weather  balloon  and  by  using  a 
second  balloon  to  host  a  passive  target  array.  Both  balloons  would  operate  in  tandem 
controlled  ascent.  Currently  C'l  profiles  are  measured  using  thermal  instrumentation. 
Using  thermal  instrumentation  restricts  operations  to  nighttime,  when  solar  loading 
on  the  instruments  is  minimized. 

The  tilt-covariance  method  also  has  several  advantages  over  the  tri-scintillilation 
method  recently  published  by  Vetelino.  et  al  [30].  The  scintillation-based  approach  is 
range-limited,  making  long-range  field  tests  cost-prohibitive  and  short-range  field  tests 
impossible  for  lack  of  detectable  scintillation.  Because  the  tilt-covariance  method  uses 
passive  extended  sources,  high-cost  sources  are  not  needed.  Optionally,  light-emitting 
diodes  (LEDs)  can  be  used  for  tests  conducted  after  twilight.  Secondly,  the  new 
approach  uses  34  nonlinear  equations  to  solve  for  only  4  unknown  parameters,  while 
the  scintillation  method  uses  only  three  equations  to  solve  for  three  unknowns.  The 
redundancy  helps  stabilize  turbulence  parameter  estimates  in  the  presence  of  noise. 

4-2  Future  Work 

As  a  natural  extension  to  this  work,  the  author  plans  to  validate  the  estimates 
presented  in  Chapter  3  by  way  of  3D  SCR  simulation  of  optical  turbulence.  Analysis 
will  include  the  effects  of  scintillation  on  the  performance  of  turbulence  parameter 
estimation.  Simulated  data  will  be  processed  through  two  concurrent  pipelines:  one 
which  ignores  scintillation,  and  one  to  include  scintillation.  Results  will  be  compared 
to  those  predicted  by  a  widely  accepted  technique  published  by  Whiteley  in  2002  [32]. 
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Additional  work  may  lead  to  predictive  modeling  of  atmospheric  seeing  conditions  as 
a  function  of  time,  date,  location,  and  weather  parameters.  Of  particular  interest  to 
astronomers  is  turbulence  parameter  profile  estimation.  Certain  star  imaging  algo¬ 
rithms  which  depend  on  crude  profile  estimates  of  non-zero  inner  scale  and  finite  outer 
scale  may  be  able  to  use  estimates  of  these  parameters  obtained  from  the  boundary 
layer.  One  such  application  is  proposed  in  [3]. 
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Appendix  A.  Approximating  Waves  in  Random  Media, 

If  we  are  to  prescribe  the  use  of  Rayleigh-Sommerfeld  (R-S)  diffraction  theory,  we 
must  first  validate  the  assumptions  for  the  case  of  diffraction  through  a  turbulent 
medium1 .  R-S  theory  was  developed  specifically  for  propagation  of  a  monochromatic 
wave  through  a  homogenous  and  isotropic  medium.  In  the  course  of  this  development, 
as  presented  by  Goodman  in  [11],  a  scalar  wave  equation  is  derived  from  Maxwell’s 
equations,  under  the  assumption  of  homogenous  permittivity,  .  By  the  relationship, 
n  =  yfc,  we  see  that  this  assumption  is  immediately  violated  because  n  is  not  ho¬ 
mogenous  in  turbulence.  Fortunately,  this  violation  is  minor  because  is  manifested  in 
a  negligible  depolarization  term  resulting  from  Earth’s  mildly  turbulent  atmosphere. 
Thus,  the  usual  scalar  wave  equation  can  be  made  to  approximate  diffraction  through 
a  turbulent  medium  [12]: 

=  0  (A.1) 

A.l  Violation  of  Helmholtz  Equation 

The  next  step  in  Goodman’s  development  of  the  R-S  diffraction  theory  necessi¬ 
tates  that  a  diverging  spherical  wave  satisfy  the  following  Helmholtz  equation  for  all 
time  [11]: 

(V2  +  [k(p)}2)9(p,t)  =  0  (A.2) 

In  this  expression,  the  propagation  constant,  k,  is  a  random  field  defined  by  the 
relationship,  k  =  and  g  is  the  candidate  Green’s  function  describing  a  diverging 
spherical  wave.  Since  the  Helmholtz  equation  is  an  entirely  spatial  relationship,  we 
define  g(p,  to)  for  some  particular  t  =  to  to  have  a  phase  offset  of  ipo'- 

g{p,to)  =  |p  -  Pol'1  exp{j[/c(p)|p  -  Pol  -  ¥*>]}.  (A.3) 

'The  medium  is  also  assumed  to  be  nonmagnetic,  dielectric,  nondispersive,  and  linear.  [11] 
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It  is  now  a  trivial  matter  to  show  that  a  diverging  spherical  wave  does  not  satisfy 
the  Helmholtz  equation  in  a  turbulent  medium  represented  by  the  random  field,  k(p). 
We  begin  by  defining  the  Laplacian  operator,  V2,  in  spherical  coordinates  [31]: 

2  A  d2  2d  1  d 2  cos(cf>)  d  Id2  .  . 

dr2  rdr  r2  sin2  (tj>)  dO2  r2  sin(0)  d<t>  r2  d<j>2' 

Next,  we  allow  k  to  be  spatially  variant  and  use  the  definition  from  (A. 4)  in  applying 
(A.'2)  to  (A. 3).  Expanding  the  left-hand  side  of  (A.'2)  yields  the  following  expression 
for  the  error  in  the  electromagnetic  field,  expressed  in  spherical  coordinates: 


Aerr  =  9 


n  ,dk  (  dk\  d2k  .dk 

2rk  dr  \dr  )  dr2  +  J  dr 


+  r  1  sin  ( <t> ) 


dk 


.d2k 


dO  )  '  J  dO2 


.  ,.dk 

jcot^)- 


fdk\ 2  d2k\ 
\d4>)  +Jd(j>2)  ’ 

(A.5) 


The  spatial  dependencies  on  g  and  k  have  been  dropped  for  convenience,  and  the 
chosen  origin  is  p0  so  that  r  —  \p  —  po  |  •  Ignoring  this  electromagnetic  field  error  is 
equivalent  to  asserting  that  Vfc  «  0. 

Arguably,  no  possible  assignment  of  g  may  exist  to  satisfy  (A. 2),  but  such  a  proof 
is  beyond  the  scope  of  this  work.  Instead  of  meticulously  analyzing  the  existence  of 
a  perfect  solution,  we  will  simply  assume  any  linear  combination  of  spherical  waves 
approximately  satisfies  the  Helmholtz  equation  in  turbulence.  This  is  nearly  equivalent 
to  the  approximation  made  by  Goodman  when  by  discarding  the  depolarization  term 
in  arriving  at  (A.l)  [11]. 
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A. 2  Turbulence- Modified  Diffraction  Formula 

The  free-space  R-S  diffraction  formula  is  given  by  [11] 

OO 

u(Po)  =  jy  IJ  exp{jk0r}  cos  (-2,  r)  dpx,  (A.6) 

—  OO 

where  u  is  the  scalar  field,  frozen  in  time,  p  =  ( x ,  y,  0)  refers  to  a  point  in  the  plane  of 
the  aperture  described  by  the  binary  window  function,  IT' (p±_),  where  px  =  (p  x.  p  y)- 
Po  =  (x0,  j/0)  zo)  is  second  point  in  free  space  described  by  the  vector  r  =  p0  —  p.  2  is 
the  unit  vector  normal  to  the  plane  of  the  aperture  pointing  toward  po,  and  ko  is 
the  free-space  propagation  constant.  Figure  A.l  illustrates  the  geometry  for  (A.6). 
The  derivation  of  the  free-space  R-S  diffraction  formula  in  (A.6)  assumes  that  k  is 


(a)  (b) 


Figure  A.l:  Rayleigh-Sommerfeld  vacuum  propagation.  This  figure  illustrates  the 
geometry  for  (A.6). 

constant  by  setting  the  gradient  of  k  to  zero. 

We  can  modify  this  formula  to  account  for  phase  distortions  associated  with 
a  turbulent  dielectric.  We  do  this  by  assuming  that  E[k]  =  k0  and  making  the 
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1 

substitution  kfar  =  f  k(p  +  s(po  -  p))ds,  which  results  in  the  following  turbulence- 
o 

modified  R-S  diffraction  formula: 


exp  \jj  k(p  +  s{p0  -  p))ds  1  dpi,  (A. 7) 


where  we  have  also  borrowed  a  portion  of  the  Fresnel  approximation  by  substituting 
r~l  cos  (—2,  r)  =  Zq1.  It  should  be  noted  that  this  approximation  is  only  valid  for 
optical  paths  short  enough  where  scintillation  is  negligible. 
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on  the  error  associated  with  SCR.  li  is  shown  empirically  that  the  theoretical  bound  may  be  used  to  predict  relative  algorithm 
performance.  Therefore,  the  covariance  delineation  theorem  is  a  powerful  tool  for  optimizing  SCR  algorithms.  For  the  third 
contribution,  we  introduce  a  new  method  for  passively  estimating  optical  turbulence  parameters,  and  demonstrate  the  method  using 
experimental  data. 
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